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ABSTRACT 


The  Rachford-Rice  objective  function  for  flash  calculations  exhibits  a  nearly 
flat  slope  across  the  two-phase  region  and  sharp  discontinuities  near  the  dewpoint. 
These  features  make  iterative  solution  procedures  sensitive  to  the  initial  estimate  of 
the  root  and  prone  to  spurious  values  if  a  correction  step  throws  the  algorithm  out¬ 
side  the  two-phase  region  or  near  the  phase  boundary. 

This  work  centers  on  the  recasting  of  the  Rachford-Rice  objective  function  into 
a  polynomial  function  of  the  vapor  fraction,  a.  The  degree  of  this  polynomial  is 
one  less  than  the  number  of  components  in  the  system  and  its  coefficients  can  be 
calculated  from  the  feed  composition  and  the  equilibrium  ratios.  A  recursive 
expression  is  developed  that  involves  symmetric  functions  and  can  be  easily  pro¬ 
grammed  on  a  computer  or  scientific  calculator. 

The  principal  advantage  of  this  new  form  of  the  objective  function  is  that  the 
theory  of  polynomials  is  well-developed.  The  location  of  their  zeros  can  be 
predicted  with  confidence  by  techniques  based  on  sound  mathematical  principles, 
such  as  the  Fourier-Budan  theorem.  The  a-polynomial  is  well-behaved  over  the 
two-phase  region  and  its  root  can  be  quickly  located  by  a  hybrid  method  of 
interval-halving  technique  and  Newton-Raphson  procedure.  The  validity  of  the  new 
objective  function  and  its  automatic  coefficient-generating  algorithm  are  tested  using 
several  multicomponent  systems  for  which  experimental  data  are  available. 

Results  of  these  tests  prove  conclusively  the  validity  of  the  generalized  polyno¬ 
mial  objective  function.  The  versatility  of  this  form  of  the  flash  objective  function, 
compared  with  the  original  Rachford-Rice  version,  is  demonstrated.  Another 


IV 


potential  advantage  of  the  polynomial  form  is  its  ability  to  handle  dilute  systems  in 
which  some  components  are  present  but  in  very  low  concentrations.  It  also  prom¬ 
ises  possible  usage  as  a  means  of  developing  appropriate  lumping  schemes. 
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=  vapor  fraction  [Equation  (1.1)] 

=  isothermal  compressibility  (Section  2.3.2) 

=  binary  interaction  coefficient  [Equation  (2.16)] 

=  bulk  modulus  [Equation  (2.23)] 

=  function  to  describe  a(T)  away  from  critical  point  [Equation  (2.11)] 
=  function  of  acentric  factor  in  PREOS  [Equation  (2.13)] 

=  function  of  acentric  factor  in  PRSV  EOS  [Equation  (2.26)] 

=  parameter  in  PRSV  EOS  [Equation  (2.27)] 

=  coefficient  of  polynomial  [Equation  (5.2)] 

=  multiplicity  of  a  part  in  a  partition  [Equation  (4.1)] 

=  EOS  variable  [Equations  (2.24),  (2.25)]] 

=  Pitzer  acentric  factor  for  the  i  -th  component  [Equation  (2.5)] 


B  =  boiling  point  (Equation  (2.1)] 

c  =  critical  property 

i  J  Jc,ljn  =  individual  components  of  the  fluid  system 

ij  =  interaction  between  component  i  and  component  j  of  the  fluid  system 
k  =  convergence  pressure  (Section  2.2) 


xn 


P  -  constant  pressure  [Equation  (4.29)] 

R  =  reduced  property 

RA  -  Rackett  compressibility  factor  [Equation  (2.22)] 

r  =  order  of  symmetric  function  [Equation  (3.41)] 

T  =  constant  temperature  [Section  2.3.2,  Equation  (4.29)] 


Superscripts 


I  =  fluid  state  [Equations  (2.24),  (2.25)] 
i  =  component  of  fluid  system  [Equations  (2.24),  (2.25)] 
(n)  =  order  of  derivative 

'  =  first  derivative 


=  second  derivative 
=  third  derivative 


Abbreviations 

CPU  =  computer  central  processing  unit 

EOS  =  equation  of  state 

°F  =  degree  Fahrenheit 

LHS  =  left-hand  side  (of  an  equation) 

PREOS  =  Peng-Robinson  equation  of  state 

PRSV  EOS  =  Peng-Robinson-Stryjek-Vera  equation  of  state 
psia  =  pounds  per  square  inch,  absolute 
Q.E.D.  =  quod  erat  demonstrandum,  which  was  to  be  proved  (Appendix  A) 
RHS  =  right-hand  side  (of  an  equation) 

SRKEOS  =  Soave-Redlich-Kwong  equation  of  state 


SSM  =  successive  substitution  method 
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Chapter  1 


DISCUSSION  OF  THE  PROBLEM 


1.1  Introduction 


Detennination  of  the  equilibrium  state  of  coexisting  liquid  and  vapor  phases, 
particularly  for  multicomponent  fluid  mixtures,  is  of  vital  interest  to  the  petroleum 
and  chemical  industries.  Many  processes  in  petroleum  production  and  refining 
involve  repetitive  flash  calculations  for  design  and  operational  purposes.  The  pri¬ 
mary  goal  of  performing  flash  calculations  is  to  determine  the  relative  amounts  and 
compositions  of  the  coexisting  phases  for  a  given  feed  composition  at  a  specified 
condition  of  temperature  and  pressure. 

This  work  is  confined  solely  to  two-phase  vapor-liquid  equilibrium  computa¬ 
tions,  although  its  results  will  no  doubt  find  application  in  multiphase  flash  prob¬ 
lems  in  the  future. 


1.2  The  Generic  Flash  Algorithm 

To  begin  the  calculation,  the  following  variables  must  be  specified:  the  system 
pressure  and  temperature,  the  molar  composition  of  the  feed  stream,  z, ,  and  an  ini- 

y. 

tial  estimate  of  the  equilibrium  ratios,  AT,  =  — .  The  process  is  assumed  to  occur 
under  isothermal  and  isobaric  conditions.  The  stages  of  the  calculation  are: 


1.  Compute  initial  estimates  of  the  equilibrium  ratios  by  one  of  the  established 
techniques  or  by  an  empirical  correlation. 
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2.  Calculate  the  phase  distribution  and  compositions  corresponding  to  the  given 
A(-values.  This  involves  the  iterative  solution  of  the  following  objective  func¬ 
tion,  developed  by  Rachford  and  Rice  (1952): 


N 

I 


i=l 


Zj  (Kj-l) 

1  +  a  (  Ki  -  1  ) 


=  0 


where  a  is  the  vapor  molar  fraction. 


(U) 


3.  Use  an  equation  of  state  (EOS)  to  calculate  the  component  fugacities  in  each 
phase  and  check  for  equality. 

4.  If  equality  is  not  achieved  (i.e.,  the  phases  are  not  in  equilibrium),  correct  the 
A'-values  on  the  basis  of  the  fugacities  and  repeat  steps  2-4. 

This  correction  is  commonly  perfoimed  using  a  successive  substitution-type 
method  or  a  second-order  Newton-type  scheme.  These  algorithms  are  well-known 
and  are  described  in  several  papers  [e.g.,  Risnes  et  al.  (1981);  Michelsen  (1982); 
Boston  and  Britt  (1978)]. 

Successful  implementation  of  the  generic  flash  algorithm  described  above 
requires  three  principal  elements.  These  are  (1)  a  general  estimate  of  the  set  of 
equilibrium  ratios  to  start  the  procedure;  (2)  a  good  equation  of  state  to  improve  AT, ; 
and  (3)  a  robust  objective  function  that  guarantees  convergence  to  a  single  value  of 
a.  A  poor  first  guess  of  AT-values  may  produce  a  phase  split  that  is  physically 
impossible  under  the  prevailing  pressure  and  temperature.  Satisfactory  methods  are 
available  for  generating  these  values.  Furthermore,  existing  equations  of  state  do  a 
fairly  good  job  of  predicting  phase  properties,  and  other  efforts  continue  along  this 
line. 
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One  area  that  has  not  enjoyed  equal  amounts  of  attention  for  a  long  time  is  the 
form  of  the  objective  function.  Invariably,  the  Rachford-Rice  objective  function 
[Equation  (1.1)]  is  most  often  used.  Recent  investigations  (Warren,  1991)  have 
shown  that  this  objective  function  does  exhibit  some  strange  behavior  which  may 
affect  its  ability  to  generate  good  results  for  some  conditions. 

1.3  Objectives  of  the  Investigation 

This  study  is  aimed  at  evolving  a  generalization  of  the  new  polynomial  form 
of  the  Rachford-Rice  objective  function  developed  by  Warren  (1991).  The  resulting 
generalized  polynomial  function  of  the  vapor  fraction,  a,  should  be  capable  of  han¬ 
dling  an  ^-component  mixture.  The  coefficients  of  the  generalized  polynomial 
should  depend  on  only  two  variables,  the  molar  composition  of  the  feed  stream  and 
the  equilibrium  ratios,  and  should  be  easy  to  obtain,  either  analytically  or  numeri¬ 
cally.  Appropriate  algorithms  are  to  be  developed  for  this  purpose. 

The  principal  advantage  of  a  polynomial  form  of  the  flash-calculation  objective 
function  is  that  the  theory  of  polynomials  is  well-developed  and  semi-analytical 
solution  techniques  exist  for  equations  up  to  fifth-order  (Zaguskin,  1961).  For 
higher-order  polynomials,  the  Newton-Raphson  iterative  method  usually  provides  a 
fast  and  accurate  determination  of  the  roots. 

Determination  of  all  the  zeros  of  this  polynomial  is  unnecessary  since  the  phy¬ 
sics  of  the  problem  demands  that  only  the  zeros  on  the  bounded  interval  [0,1]  are  of 
practical  interest.  Furthermore,  the  physics  also  suggests  that  only  one  zero  (or 
value  of  a)  exists  on  this  interval,  which  represents  the  two-phase  vapor-liquid 
region.  It  can  be  shown  mathematically  that  this  is  indeed  the  case  for  well-defined 
systems,  as  will  be  demonstrated  in  52. 


Chapter  2 


LITERATURE  REVIEW 

A  survey  of  the  pertinent  literature  reveals  that  apparently  only  one  other 
worker,  Warren  (1991),  has  studied  the  particular  aspect  of  flash  calculations  tar¬ 
geted  in  this  research.  A  comprehensive  review  of  the  literature  pertaining  to  the 
use  of  cubic  equations  of  state  in  flash  calculations  was  conducted  in  order  to  pro¬ 
vide  a  reference  point  for  the  testing  of  the  polynomial  objective  function. 

This  review  is  sub-divided  into  three  sections:  flash  calculation  algorithms; 
equilibrium  ratios;  and  cubic  equations  of  state.  Particular  emphasis  is  laid  on  the 
Peng-Robinson  equation  of  state. 

2.1  Vapor-Liquid  Equilibrium  Flash  Calculations 

This  discussion  will  be  confined  to  two-phase  vapor-liquid  equilibria.  The 
work  to  date  concentrates  on  developing  robust  algorithms  with  rapid  convergence 
rates.  Robustness  implies  the  ability  to  continue  the  calculations  after  recovering 
from  a  spurious  value  of  the  vapor  fraction  computed  in  the  neighborhood  of  the 
critical  point  or  at  the  phase  boundaries.  Abhvani  and  Beaumont  (1987)  present  an 
excellent  review  of  EOS-based  flash  algorithms.  They  divide  the  papers  into  two 
categories  according  to  solution  method,  those  using  some  variant  of  the  successive 
substitution  method  (SSM)  or  those  employing  a  second-order  Newton-type  method. 

The  SSM  technique  is  the  traditional  solution  algorithm,  but  it  exhibits  a  poor 
rate  of  convergence  and  has  stability  problems  close  to  saturation  points  and  in  the 
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critical  region.  Risnes  et  al.  (1981),  Michelsen  (1982),  and  Mehra  et  al.  (1983) 
made  attempts  at  acceleration  and  stabilization  of  this  method. 

Similarly,  many  workers  have  proposed  various  forms  of  second-order  Newton 
procedures  to  avoid  the  slow  rate  of  convergence  of  the  SSM,  such  as  Boston  and 
Britt  (1978),  Fussell  and  Yanosik  (1978),  Asselineau  et  al.  (1979),  Fussell  (1979), 
Baker  and  Luks  (1980),  and  Varotsis  et  al.  (1981).  Others  advocate  a  combination 
of  successive-substitution  and  Newton  methods;  the  former  is  used  to  provide  good 
initial  values  to  the  rapidly  converging  latter.  Informative  studies  include  Mon 
(1980,  1983),  Mehra  et  al.  (1982),  Michelsen  (1982),  Nghiem  et  al.  (1983),  and 
Abhvani  and  Beaumont  (1987). 

Benmekki  (1984)  developed  a  general  algorithm  for  flash  calculations  that  can 
utilize  any  cubic  equation  of  state  and  features  a  specified  calculational  path  for 
computing  the  phase  boundaries.  This  is  an  attempt  to  ensure  that  bubblepoint  and 
dewpoint  computations  originate  from  within  the  two-phase  region,  thus  guarantee¬ 
ing  meaningful  values  of  the  equilibrium  ratios. 

Warren  (1991)  made  a  radical  departure  from  previous  efforts  at  enhancing 
flash  calculation  algorithms  when  he  formulated  an  explicit  linear  equation  for  the 
vapor  fraction  of  a  binary  system.  He  successfully  extended  this  to  a  quadratic 
equation  for  a  ternary  system  and  a  cubic  equation  for  a  quaternary  mixture.  The 
success  achieved  by  Warren  and  the  possibility  of  the  existence  of  a  generalized 
polynomial  expression  for  the  vapor  fraction  in  a  two-phase,  N-component  fluid  sys¬ 
tem  motivated  the  current  work. 
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2.2  Vapor-Liquid  Equilibrium  Ratios 

The  use  of  initial  equilibrium  ratios  close  to  the  final  values  for  a  multicom¬ 
ponent  fluid  is  crucial  to  the  rapid  convergence  of  any  flash  calculation.  Experi¬ 
mental  values  are  preferred  because  the  prediction  of  AT,  for  a  particular  fluid  at 
various  combinations  of  temperature,  pressure  and  composition  requires  lengthy  cal¬ 
culations.  Therefore,  predictive  methods  for  AT-values  are  a  limiting  factor  in  the 
speed  and  robustness  of  any  flash  calculation  algorithm. 

The  expression  "equilibrium  constant"  was  coined  by  Souders  et  al.  (1932) 
and  was  defined  as  the  ratio  of  the  vapor  mole  fraction  to  that  of  the  liquid.  The 
basis  for  most  predictive  methods  had  its  genesis  when  Cox  (1923)  observed  that 
the  lines  on  a  semilogarithmic  plot  of  vapor  pressure  against  temperature  appeared 
to  converge  to  a  single  pressure.  Katz  arid  Hachmuth  (1937)  demonstrated  an 
analogous  behavior  for  equilibrium  constants;  they  converged  to  unity  at  a  fluid 
mixture’s  critical  pressure. 

White  and  Brown  (1942)  attempted  to  develop  a  correlation  for  Ar-values  based 
on  this  "convergence"  pressure.  Hanson  and  Brown  (1945)  used  experimental  data 
to  correlate  the  convergence  pressure  (/>*)  at  one  temperature  as  a  function  of  the 
molal  average  boiling  point  of  the  equilibrium  vapor  and  liquid.  They  also  showed 
that  the  convergence  pressure  concept  could  be  extended  from  binary  to  multicom¬ 
ponent  systems. 

Hadden  (1948,  1953)  produced  nomographs  for  equilibrium  constants  of  pure 
components  as  functions  of  temperature  and  pressure,  and  incorporated  convergence 
and  vapor  pressures  into  nomographs  for  mixtures.  He  demonstrated  that  mixture 
convergence  pressure  is  a  function  of  the  operating  temperature  and  of  the  liquid- 
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phase  composition  exclusive  of  the  lightest  component  concentration.  This  compo¬ 
sition  dependence  led  Muskat  (1949)  to  propose  the  use  of  the  term  "equilibrium 
ratio"  in  place  of  "equilibrium  constant."  Edmister  (1949)  presented  a  graph  involv¬ 
ing  the  ratio  of  differences  between  the  convergence  and  critical  pressures  and  the 
ratio  of  differences  between  the  system  and  critical  temperatures. 

Winn  (1952)  developed  nomographs  based  on  Hadden’s  (1948)  results  that 
allow  the  determination  of  ^-values  at  a  convergence  pressure  of  5000  psia.  For 
systems  with  pk  *  5000,  he  provides  a  translation  to  find  the  value  of  Kt  at  other 
"apparent"  convergence  pressures.  The  methods  proposed  by  these  three  authors 
require  charts  and  do  not  lend  themselves  to  computer  calculations. 

Hoffmann  et  al.  (1953)  attempted  to  extend  Cox’s  (1923)  vapor  pressure  graph 
for  the  purpose  of  determining  equilibrium  ratios  by  plotting  log  Kp  against  the 
component  characterization  factor  F,  where 


Kp  -  product  of  equilibrium  ratio  and  pressure 


b  -  constant  required  to  translate  the  vapor  pressure  curve 
for  a  hydrocarbon  onto  the  straight  line  of  the  Cox  chan 
Tb  =  hydrocarbon  boiling  point 
T  -  system  temperature 


Brinkman  and  Sicking  (1960)  presented  an  iterative  method  for  finding  the  conver¬ 
gence  pressure  based  on  the  slope,  sP ,  of  the  plot  mentioned  in  Hoffmann  et  al. 
(1953).  Then,  the  equilibrium  ratio  could  be  determined  as 


*  =  ^ i- 

P 


(2.2) 


8 


Standing  (1979)  observed  that  the  composition  dependence  of  the  equilibrium 
ratio  is  negligible  at  pressures  below  1000  psia.  He  proceeded  to  combine  the  work 
of  Hoffmann  et  al.  (195?)  and  Brinkman  and  Sicking  (1960)  to  develop  a  correla¬ 
tion  for  K-values  for  the  crude  oils  studied  by  Katz  and  Hachmuth  (1937): 


tf  =  llO<a+cF>  (2.3) 

P 

where  a  and  c  are  the  intercept  and  slope  (respectively)  of  log  Kp  vs.  F  plots  of  the 
abovementioned  oils.  Both  a  and  c  are  expressed  as  functions  of  pressure.  He  also 
presented  equations  for  the  heavy  fraction  and  the  common  reservoir  gases  N2,  C02 
and  H2S  (referred  to  as  permanent  gases). 

Wilson  (1969)  published  a  A'- value  equation  that  currently  enjoys  widespread 
use  in  flash  calculations: 


where 

B  =  5.37(1  +to,)(l--^-) 

‘  Ri 

pRi  =  reduced  pressure  of  the  i  -th  component 
TRi  =  reduced  temperature  of  the  i  — th  component 
Oty  =  Pitzer  acentric  factor  of  the  i  -th  component 


(2.4) 

(2.5) 


Wilson’s  equation  fails  to  predict  accurate  equilibrium  ratios  for  most  fluids  above 
pressures  of  500  psia,  as  illustrated  by  Warren  (1991).  Whitson  and  Torp  (1981) 
attempted  to  correct  this  problem  by  re-introducing  the  system  convergence  pressure 
to  the  Wilson  equation: 


= 


Pci 

Pk 

PRi 

(2.6) 
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where 


A  =  1  - 


JL 

Pk 


14.7 


14.7 


0.6 


(2.7) 


pci  =  critical  pressure  of  the  i  -th  component 

Risnes  and  Dalen  (1984)  took  an  approach  based  on  the  equation  of  state  used 
in  the  flash  calculations.  Their  basic  idea  was  to  assume  the  mixture  or  feed  to  be 
liquid  and  then  evaporate  up  to  one-half  of  the  system  to  form  a  gas  phase  by  use 
of  the  fugacities.  The  initial  ^-values  then  could  be  calculated  from  the  resulting 
phases.  This  method  is  reported  to  perform  well  near  the  critical  point  and  along 
the  bubblepoint  line  but  often  fails  along  the  dewpoint  curve. 

Reportedly,  the  most  accurate  lvalue  predictor  is  that  proposed  by  Varotsis 
(1989).  He  used  over  1000  experimental  equilibrium  ratios  to  construct  an  X-Y  plot 
such  that  each  reservoir  fluid’s  position  on  the  "map"  is  determined  by  its  coordi¬ 
nates  X  and  Y.  These  coordinates  are  described  by  a  polynomial  fitted  to  the 
apparent  pressure  mentioned  in  Winn  (1952).  He  proposes  an  equation  for  the  con¬ 
vergence  pressure  based  on  the  mole  fraction  of  methane  and  nitrogen  in  the  fluid. 

Each  pure  hydrocarbon  component  is  represented  on  the  map  by  its  own  set  of 
coordinates  {Xt ,  ),  which  are  calculated  as  functions  of  the  component  acentric 

factor.  Specific  values  are  given  for  the  permanent  gases  and  correlations  based  on 
molecular  weight  are  specified  for  the  heptane-plus  fraction. 

Finally,  the  straight  line  that  joins  the  pressure  and  temperature  coordinates 
(X,  T)  of  the  fluid  with  the  position  of  each  component  on  the  map  (Xiy  T,  )  inter¬ 
sects  the  AT-value  axis  at  a  point  that  corresponds  to  the  equilibrium  value  of  the 
selected  constituent.  Varotsis  (1989)  presents  tables  for  three  different  crude  oils 
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and  gas  condensates  at  varying  temperatures  and  pressures  that  display  AT-values 
remarkably  close  to  experimental  values.  They  are  an  order-of- magnitude  improve¬ 
ment  over  those  predicted  by  the  equations  of  Wilson  (1969)  or  Whitson  and  Toip 
(1981). 

The  method  of  Varotsis  (1989)  was  attempted  in  the  current  work.  His  K- 
value  predictor  was  formulated  using  data  for  crude  oils  and  gas  condensates  con¬ 
taining  the  Ci  -  C  6  alkane  series,  the  heptane-plus  pseudocomponent  and  the  per¬ 
manent  gases.  It  will  not  properly  describe  systems  (such  as  the  methane-ethane- 
propane  ternary)  containing  fewer  components  than  these  "typical"  reservoir  fluids. 
For  lack  of  a  suitable  replacement  expression  for  pk  ,  Wilson’s  equation  is  used  in 
the  cun-ent  work. 

2.3  Cubic  Equations  of  State  (EOS) 

The  equation  of  state  (EOS)  is  the  heart  of  a  modem  flash  calculation  algo¬ 
rithm.  Ideally,  it  should  be  able  to  accurately  represent  the  thermodynamic  proper¬ 
ties  of  the  fluid  of  interest  over  the  complete  range  of  operating  pressures  and  tem¬ 
peratures.  Since  engineering  applications  rarely  focus  on  an  isolated  chemical 
species,  the  EOS  should  incorporate  mixing  rules  that  allow  it  to  extend  its  predic¬ 
tive  capabilities  to  the  behavior  of  multicomponent  fluids.  Its  component-specific 
descriptive  parameters  should  be  readily  calculable  from  well-known  properties, 
such  as  critical  temperature  and  pressure,  molecular  weight  and  acentric  factor. 
Finally,  the  associated  computations  should  not  consume  excessive  computer  time, 
especially  if  the  equation  of  state  is  to  be  used  for  repetitive  calculations. 

The  engineer  is  faced  with  the  choice  of  using  a  complex  EOS  exhibiting  a 
high  degree  of  non-linearity  and  many  adjustable  parameters,  or  a  cubic  EOS  which 
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possesses  an  analytical  solution  and  requires  the  estimation  of  two  or  three  parame¬ 
ters.  Mathias  and  Benson  (1986)  presented  a  comparison  of  average  central¬ 
processing-unit  (CPU)  times  required  by  three  cubic  EOS  and  by  three  complicated 
EOS  to  compute  fugacity  coefficients  and  enthalpy  departures.  They  asserted  that 
the  time  required  for  any  of  the  candidate  equations  to  calculate  the  density  root  (or 
compressibility  factor  root)  is  negligible  compared  to  that  involved  in  executing  the 
various  mixing  rules.  In  fact,  for  systems  containing  more  than  about  six  com¬ 
ponents,  the  cubic  EOS  become  more  computationally  burdensome  than  the  compli¬ 
cated  EOS  simply  because  of  the  cross  terms  inherent  to  the  cubic  EOS  mixing 
rules. 

Engineers  frequently  use  cubic  EOS  because  they  work  well  over  the  range  of 
most  industrial  operating  conditions  and  are  easily  programmed  for  solution  on  a 
computer.  The  two  cubic  EOS  which  have  gained  the  widest  acceptance  are 
Soave’s  modifications  of  the  Redlich-Kwong  (1949)  equation  of  state  (SRKEOS) 
(Soave,  1972)  and  that  presented  by  Peng  and  Robinson  (1976b)  (PREOS).  The 
PREOS  and  suggested  improvements  are  examined  in  this  work  for  possible  use  in 
flash  calculations  because  of  the  author’s  familiarity  with  this  EOS. 

23.1  Development  of  the  Peng-Robinson  EOS 

Upon  the  success  of  the  SRKEOS,  Peng  et  al.  (1975)  undertook  a  further  study 
to  formulate  a  cubic  equation  of  state  with  an  improved  capability  to  predict  liquid 
densities  and  other  fluid  properties,  particularly  in  the  vicinity  of  the  critical  region. 

This  study  resulted  in  a  further  modification  of  the  attractive  pressure  term  of 
the  classical  equation  of  state  proposed  by  van  der  Waals  (1873).  The  result  was 
the  EOS  presented  by  Peng  and  Robinson  (1976b): 
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_  RT _ g(T) 

P  v-b  v(v+b)  +  b(v-b) 

The  values  of  the  parameters  are  obtained  from 


b  =  0.07780  — -  (2.9) 

Pc 

Rip  2 

a  (Tc )  =  0.45724  - —  (2. 10) 

Pc 

a<T)  =  a(Tc)T\(TR,<fi)  (2.11) 

r\A  =  1  +  k(1  -TRA)  (2.12) 

k  =  0.3746  +  1.48500)  -  0.1644o)2  +  0.01667a)3  (2.13) 


Equation  (2.12)  has  the  same  form  as  that  used  by  Soave  (1972)  but  k  was  obtained 
by  fitting  a  larger  range  of  vapor  pressure  data  as  a  function  of  the  reduced  tem¬ 
perature  and  the  acentric  factor  (Pitzer  et  al.t  1955)  of  each  component. 

In  order  to  use  the  equation  for  systems  containing  more  than  one  component, 
the  following  mixing  rules  are  presented: 


a  =  XX  xixjaij  (2.14) 

*  j 

b^Zxibt  (2.15) 

i 

where 

Oij  =d  -SfjfrHaj*  (2.16) 

Equations  (2.14)  and  (2.15)  are  consequences  of  the  mixing  rule  proposed  by  Kay 
(1936),  while  Equation  (2.16)  was  developed  by  Zudkevitch  and  Joffe  (1970).  The 
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experimentally  determined  binary  interaction  coefficient,  btJ ,  characterizes  the 
binary  formed  by  component  i  and  component  j.  The  importance  of  btJ  in  accu¬ 
rately  reproducing  P-V-T  data  was  discussed  by  Peng  and  Robinson  (1976b)  and  by 
Robinson  et  al.  (1985). 

The  PREOS  can  be  written  in  the  form  of  a  cubic  equation  in  the  compressibility 
factor: 

Z3  -  (1  -B)Z2  +  (A  -  3B2  -2B)Z  -  (AB  -  B2  -  B3)  =  0  (2.17) 

where 


A  = 


_S2_ 

r2t2 


(2.18) 


B  =  -|p  (2.19) 

z  =  (2.20) 

2.3.2  Selection  of  the  Proper  Root  in  Cubic  EOS 

Equation  (2.17)  yields  one  or  three  roots  depending  upon  the  number  of  phases 
in  the  system.  The  authors  stated  that,  in  the  two-phase  region,  the  largest  root  is 
for  the  compressibility  factor  of  the  vapor  while  the  smallest  positive  root 
corresponds  to  that  of  the  liquid. 

Lawal  (1987),  however,  asserted  that  this  criterion  was  insufficient  to  select  the 


proper  root.  He  proved  that,  in  the  event  of  multiple  real  roots,  the  smallest  of  the 
positive  roots  larger  than  or  equal  to  B  must  be  chosen  for  the  compressibility  of 
the  liquid. 
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Asselineau  et  al.  (1979)  compared  the  calculated  volume  to  the  pseudo-critical 
volume  to  assign  the  root  to  the  proper  phase,  under  specific  conditions.  Poling  et 
al.  (1981)  examined  the  order  of  magnitude  of  the  isothermal  compressibility, 
P  =  —0v  /dp  )f  Iv ,  to  ascertain  the  presence  of  the  liquid  or  vapor  phase.  Gosset  et 
al.  (1986)  offered  two  discriminants,  one  based  on  the  Cardan  criterion  for  the 
number  of  real  roots  for  a  cubic  equation  and  a  heuristic  approach  similar  to  that  of 
Asselineau  et  al.  (1979). 

2  J.3  Modifications  to  the  Peng-Robinson  EOS 

Numerous  attempts  have  been  made  to  correct  for  the  deficiencies  inherent  in  a 
cubic  equation  of  state  by  introducing  additional  parameters  into  the  PREOS. 
These  changes  improve  some  aspect  of  the  EOS’s  performance  (usually  liquid  den¬ 
sity  predictions)  but  at  the  cost  of  increased  complexity  and  the  requirement  for 
tables  or  correlations  to  determine  the  additional  parameter(s)  for  each  fluid  com¬ 
ponent.  This  review  will  touch  on  a  limited  number  of  these  studies. 

2.3.3. 1  Volume  Corrections 

The  modification  of  the  SRKEOS  proposed  by  Peneloux  et  al.  (1982)  also 
formed  the  basis  for  two  other  studies  concerned  with  the  PREOS.  These  authors 
suggested  that  the  use  of  a  "pseudo  volume"  defined  by 

v  =  v  +  2  cixi  (2.21) 

i 

could  be  used  to  effect  a  translation  along  the  volume  axis,  leaving  unchanged  the 
predicted  equilibrium  conditions.  They  chose  c  so  that  correct  saturated  liquid 
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densities  were  exactly  reproduced  at  the  reduced  temperature  TR  =0.7.  They 
rejected  the  acentric  factor  as  a  correlating  parameter  in  favor  of  the  Rackett 
compressibility  factor,  Zra  .  developed  by  Spencer  and  Danner  (1972): 


RT 

c  =  0.40768 — —(0.29441  -  )  (2.22) 

Pc 

Their  third  parameter  did  improve  predictions  of  saturated  liquid  densities. 

Almost  simultaneously,  Jhaveri  and  Youngren  (1988)  and  Mathias  et  al. 
(1989)  presented  three-parameter  modifications  of  the  PREOS  based  on  the  work  of 
Peneloux  et  al.  (1982).  The  first  authors  correlated  the  third  parameter,  c,  with 
molecular  weight.  The  second  study  retained  the  Pdneloux-Rauzy-Freze  volume 
correction  scheme  but  added  a  further  term  involving  the  bulk  modulus  to  handle 
the  critical  region.  The  bulk  modulus  is  dimensionless  and  is  defined  as: 


5 


(2.23) 


From  an  examination  of  the  graphs  accompanying  both  publications,  the  work  of 
Mathias  et  al.  (1989)  seems  to  produce  results  closer  to  the  experimental  values  for 
saturated  volumes  and  densities. 


23.32  Temperature  Dependence 

Xu  and  Sandler  (1987a,b)  postulated  that  the  molar  co-volume  term,  b,  is  not 
independent  of  temperature  and  they  disputed  the  fitting  of  vapor  pressures  used  by 
Peng  and  Robinson  (1976b)  to  characterize  the  attractive  constant,  a.  They  corre¬ 
lated  vapor  pressure  and  volume  data  for  16  components  at  both  subcritical  and 
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supercritical  conditions  and  proposed  to  replace  the  numeric  coefficients  of  a  and  b 
found  in  Equation  (2.9)  and  Equation  (2.10)  with: 

(2.24) 

i=0 

and 

Vt,1  =  b'Tk  (2.25) 

«= o 

where  i  refers  to  the  species  and  I  denotes  either  subcritical  or  supercritical  condi¬ 
tions. 

Wu  and  Sandler  (1989)  generalized  the  temperature-dependent  parameters  of 
Xu  and  Sandler  (1987a,b)  by  performing  least-squares  fits  of  ya  and  \\rb  as  func¬ 
tions  of  acentric  factor  and  reduced  temperature.  They  were  able  to  accomplish  this 
task  only  for  the  n-alkane  series  because  of  insufficient  data.  For  their  intended 
application  of  the  work  (petroleum  reservoir  simulation),  they  envisioned  the  use  of 
the  fluid-specific  parameters  for  the  permanent  gases,  water  and  light  ends  and  the 
generalized  parameters  for  the  heavy  pseudocomponents. 

Stryjek  and  Vera  (1986a, b,c,d)  re-worked  Equation  (2.13)  to  obtain  a  better 
reproduction  of  vapor  pressure  data  at  low  reduced  temperatures: 

Kq  =  0.378893  +  1.4897153o>  -  0.17131848co2  +  0.0196554m3  (2.26) 

and  modified  Equation  (2.12)  by  the  introduction  of  one  compound-characteristic 
adjustable  parameter,  Kj: 

k  =  kq  +  Kj  (1  +Tr  *)(0.7  -  Tr  )  (2.27) 

Stryjek  and  Vera  (1986b)  and  Proust  and  Vera  (1989)  listed  values  of  Kj  for  over 
160  compounds  of  industrial  interest.  Stryjek  and  Vera  (1986d)  and  Wilczek-Vera 
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and  Vera  (1987)  examined  mixing  rules  of  varying  complexity  for  use  with  the 
Peng-Robinson-Stryjek-Vera  (PRSV)  EOS.  For  the  current  work,  the  PRSV  EOS 
with  the  original  PREOS  mixing  rules  (as  formulated  by  Stryjek  and  Vera,  1986b) 
is  used  and  produces  noticeably  better  results  than  did  the  PREOS. 

2.3.4  References  on  Cubic  EOS 

Abbott  (1979)  and  Martin  (1979)  presented  comprehensive  reviews  of  cubic 
equations  of  state  available  at  that  time,  and  Vidal  (1983)  and  Vera  et  al.  (1984) 
updated  the  topic.  Huron  and  Vidal  (1979)  proposed  composition-dependent  mixing 
rules  while  Mathias  and  Copeman  (1983)  discussed  mixing  rules  dependent  on 
volume.  Finally,  Peng  and  Robinson  (1976b),  Peng  and  Robinson  (1977),  Robinson 
and  Peng  (1978),  Robinson  (1979)  and  Peng  (1986)  developed  specific  applications 
of  their  EOS. 


Chapter  3 


DEVELOPMENT  OF  THE  POLYNOMIAL  FUNCTION 
FOR  SIMPLE  SYSTEMS 


This  chapter  discusses  the  work  published  by  Rachford  and  Rice  (1952)  and 
Warren  (1991)  on  performing  flash  calculations.  It  shows  the  development  of  the 
Rachford-Rice  objective  function  [Equation  (1.1)]  and  extends  Warren’s  work  as  a 
precursor  to  developing  a  generalized,  multicomponent  equation  for  the  vapor  frac¬ 
tion. 

3.1  The  Rachford-Rice  Flash  Objective  Function 

We  will  briefly  examine  the  derivation  of  the  Rachford-Rice  objective  function 
that  is  universally  used  today  in  flash  calculations.  After  plotting  its  behavior,  it 
will  become  plain  why  it  is  so  difficult  to  solve  by  iterative  techniques  such  as  the 
Newton-Raphson  method. 

3.1.1  The  Material  Balance  Development 

Flash  calculations  are  used  to  determine  the  compositions  and  quantities  of  the 
vapor  and  liquid  phases  at  equilibrium  which  result  when  an  A-component  fluid  of  a 
particular  composition  is  subjected  to  a  particular  pressure  and  temperature.  The 
composition  of  the  feed  stream,  F,  is  denoted  by  I  z,  and  it  flashes  into  L  moles  of 
liquid  with  composition  I  Jt,  ,  and  V  moles  of  vapor  with  composition  I  y, .  The 
resulting  material  balance  equations  are: 
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F  =  L  +  V 


(3.1) 


Fzi  =  Lx,  +  Vy,- 


(3.2) 


As  defined  in  Chapter  1,  the  equilibrium  ratio  is: 


(3.3) 


and,  rearranging,  one  is  left  with  the  equation: 


yi  =  *«*/ 

Substituting  Equation  (3.4)  into  Equation  (3.2)  yields: 


(3.4) 


Fz,  =  Vx;Ki  +Lx; 

Simplify  by  isolating  the  term  and  dividing  through  by  F: 


(3.5) 


Zi  =  Xi 


(3.6) 


Dividing  Equation  (3.1)  through  by  F  and  solving  for  —  yields: 

F 


(3.7) 
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Substituting  Equation  (3.7)  in  Equation  (3.6)  and  simplifying  the  equation  results  in: 


= 


1  +  j  (  Ki  -  1  ) 


(3.8) 


N 

Imposing  the  constraint  of  £  -  1  on  Equation  (3.8)  leaves: 

i=l 


1  +  j  ( K( ■-  1  ) 


(3.9) 


Rearranging: 


N 

z 


i= 1 


z» 

1  +  J  ( AT,  -  1  ) 


1  =  0 


Recalling  Equation  (3.4),  we  can  write: 


(3.10) 


yi 


1  +  jiKi-l) 


N 

Imposing  the  constraint  of  £  =  1  on  Equation  (3.11)  yields: 

i=i 


(3.11) 


N 

I 

i=l 


Z'Kt 


1  =  0 


(3.12) 
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Combining  Equation  (3.10)  and  Equation  (3.12)  leaves: 


N 

I 


«=1 


z,  (  Kj  -  1  ) 

1  +  j  (  Ki  -  1  ) 


=  0 


Defining  the  vapor  fraction,  a,  as: 


(3.13) 


(3.14) 


and  substituting  Equation  (3.14)  into  Equation  (3.13)  yield  the  Rachford-Rice  objec¬ 
tive  function: 


N 

z 


;= i 


Zj  (  -  1  ) 


(3.15) 


3.1.2  A  Graphic  Representation  of  the  Rachford-Rice  Objective  Function 


As  the  authors  noted,  their  formulation  of  the  objective  function  was  prone  to 
instability  near  the  values  of  a  that  represented  the  phase  boundaries,  namely,  0  and 
1.  They  showed  that  the  slope  of  the  function  near  these  points  may  be  quite  steep. 
It  is  this  feature  that  tends  to  throw  derivative-based  root-finding  techniques  out  of 
the  two-phase  region,  yielding  spurious  roots. 

Figure  3.1  depicts  the  behavior  of  the  objective  function  over  a  wide  range  of 
a  for  a  binary  system  of  70%  methane  and  30%  ethane  (Bloomer  et  al.,  1953). 
Figure  3.2  does  the  same  for  a  ternary  system  consisting  of  85%  methane,  10% 
ethane  and  5%  propane  (Parikh  et  al.,  1984).  Although  values  of  the  vapor  fraction 


24 


have  no  physical  meaning  outside  the  interval  [0,1],  these  graphs  serve  to  illustrate 
how  ill-behaved  the  objective  function  is.  Its  slope  is  almost  flat  as  it  traverses  the 
two-phase  region  and  it  is  plagued  by  spiky  singularities. 

This  work  will  attempt  to  develop  a  new  expression  for  a,  one  that  possesses 
reasonable  slope  over  the  desired  interval  and  has  no  discontinuities  near  the  phase 
boundaries. 

3.2  Warren’s  Explicit  Equations  for  the  Vapor  Fraction 

Warren  (1991)  expanded  the  Rachford-Rice  objective  function  into  a  polyno¬ 
mial  in  a  for  a  binary,  ternary  and  quaternary  fluid  system.  He  did  this  by  setting 
N  equal  to  2,  3  or  4,  respectively,  and  reducing  the  resulting  equations  to  their  sim¬ 
plest  polynomial  form  by  algebraic  manipulations.  To  demonstrate  the  validity  of 
his  work,  Warren  also  showed  that  the  higher-order  polynomials  would  reduce  to 
those  for  smaller  systems  when  the  appropriate  mole  fractions  and  equilibrium  con¬ 
stants  were  removed. 

We  will  assume  (as  did  Warren)  that,  under  isobaiic  and  isothermal  conditions, 
the  equilibrium  constant  does  not  change  such  that  the  quantity  (Ki  -  1 ),  which 
appears  in  the  objective  function,  can  be  represented  by  a  constant,  C, . 

We  will  reproduce  the  entire  process  here  for  a  binary  system  but  will  show 
only  the  final  result  for  a  ternary  and  quaternary  system,  since  the  algebra  can  be 
quite  tedious  and  repetitive. 


25 


3.2.1  Binary  System 


Starting  with  the  objective  function  as  defined  in  Equation  (3.15): 


N 


,Ti  1  +  «  (  K, :  -  1  ) 


(3.15) 


and  defining  C,  =  Kt  -  1 ,  Equation  (3.15)  can  be  rewritten  as: 


N 


zi  Ci 

Y  — ! — ! —  =  o 

h  1+aC, 


(3.16) 


For  a  two-component  system,  setting  N  =  2  in  Equation  (3.16)  and  expanding 
term-wise  yields: 


z\  Cl  +  z2  ^2  _  q 


1+aCj  1  +  a  C2 
Moving  the  terms  with  the  subscript  "2"  to  the  right-hand  side  of  the  equation: 


(3.17) 


z j  Cj 


z2  Cj 


1  +  aC,  1  +  a  C2 

By  multiplying  each  side  by  (1  +a  Cj)  (1  +a  Cj),  one  obtains: 


(3.18) 


(zi  Cj)  (1  +  a  C2)  -  -  (z2  C 2)  (1  +  a  Cj) 
Expanding  each  side  yields: 


(3.19) 


ZjCj  +  aziCiC-?-  )  C  -y  ~  ctz-jC^C 


1  u  1 


2 2 


■2  u  1  ^2 


(3.20) 
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Move  a  terms  to  the  left-hand  side  of  the  equation  and  all  remaining  terms  to  the 
right-hand  side,  then  recall  that,  for  a  binary  system,  zx  +  z2  =  1 : 


ctfCjCj)  —  —  (z\Ci  +  Z2Cj)  (3.21) 

Dividing  both  sides  through  by  Cj  C2  and  substituting  ( Kt  -  1)  yields  the  explicit 
form  of  the  objective  function  for  a  binary  system: 


a  =  - 


(3.22) 


3.2.2  Ternary  System 


a"  +  a 


3 

r  ■> 

3 

2*,- 

IC; 

i=i 

j*i 

3 

n 

i*i 

_ 

3  zi 

*  l-r- 
,=1  n  q 

j*> 


=  0 


3.2.3  Quaternary  System 


(3.23) 


a3  +  a2 


4  (l-rt) 
M  Ci 


+  a 


4 

i=l 


4 

I  Cj 

u* _ : 


n  ^ 

j*> 


+ 

,=1  n  ^ 

j** 


=  0 


(3.24) 


3 3  Extension  of  Warren’s  Work  to  Larger  Systems 


Warren’s  method  can  be  used  to  develop  polynomial  expressions  for  systems 
having  five,  six  and  seven  components.  It  will  be  observed  that  the  terms  of  the 
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equation  expand  in  a  regular  fashion,  thereby  suggesting  the  possibility  of  develop¬ 
ing  a  recursive  relationship  dependent  only  on  N,  the  number  of  components.  Only 
the  algebra  for  the  five-component  system  will  be  presented,  as  those  for  six-  and 
seven-component  systems  follow  the  same  procedure. 


3.3.1  Quinary  System 


We  begin  with  the  objective  function: 


5 

I 


zi  c, 


=  0 


fii  i +  «<:,• 

which,  when  expanded  for  five  components,  becomes: 


Z  i  C  j  z2^2  z3^3  z4^4  Z5^S  ~ 

+  t - z —  +  z - z —  +  — - r —  +  z - z —  —  U 


(3.16) 


1  +  aC,  l+aC2  1  +  a  C3  l+aC4  1  +  a  C5 


Multiplying  through  by  J~[  (1  +  a  C,)  yields: 

»= l 


(3.25) 


z\  Cj(  1  +  a  C 2)  (1  +  a  C3)  (1  +  a  C 4)  (1  +  a  C 5)  + 

z2  C 2  (1  +  a  Ci)  (1  +  a  C3)  (1  +  a  C 4)  (1  +  a  C5)  + 

z3  C3  (1  +  a  Cj)  (1  +  a  C2)  (1  +  a  C4)  (1  +  a  C5)  + 

z4  C4  (1  +  a  (1  +  a  C2)  (1  +  a  C3)  (1  +  a  C5)  + 

z5  C5  (1  +  a  Cj)  (1  +  a  C2)  (1  +  a  C3)  (1  +  a  C4)  =  0  (3.26) 

Expanding  each  term: 


a 4  2 j  Cj  C2  C3  C4  C5  + 


a?  z i  C ]  (C 2  C  3  C4  +  C  2  C  3  C  5  +  C2  C  4  C5  +  C  3  C 4  C5  )  + 

a^ZjCi(C2C3  +  C2C4  +  C2C5  +  C3C4  +  C3C5  +  C4C5)  + 

CX  Z 1  C  j  (C  2+C3  +  C4+C5)  +  Zj  Cj  + 
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a4  Z2C1  C  2  C  3  C4C5  + 


tt2  Z  2  C  2  1C3  +  C1C4  +  C1C5  +  C3C4  +  C3C5  +  C4C5)  + 

(X  Z  2  C  2  (C  j  +  C3  +  C4  +  C5)  +  Z2C2  + 


a4Z3Cj  C2C3C4C5  + 

tt^  Z  3  C  3  (C  j  C  2  C  4  +  C  \  C  2  C  5  +  C  l  C4C^'^C2^'4^--^)'b 
ft2  Z  3  C  3  (C  |C2  +  CjC4+CjC'5  +  C2C'4  +  C2C5  +  C4C5)  + 
az3  C 3  (Ci  +  C2  +  C4  +  C5 )  +  Z3  C3  + 


Z4C1C2C3C4C5  + 


a2  z4  C4  (C!  C2  +  C1C3  +  C1C5  +  C2C3+  C2  C5  +  C3  C5  )  + 

Ot  Z  4  C  4  (C  i  +  C2  +  C3  +  C5)  +  Z4C4  + 
ft4  Z5C1  C2C3C4C5  + 

CX^  Z  5  C  5  (C 1  C2C3  +  C1  C2C4  +  C1  C3  C4  +  C2C3C4)  + 
ft2  z5  C5  (C 1C2  +  C1C3  +  C1C4  +  C2C3  +  C2C4  +  C3C4)  + 


=  0 


(3.27) 


5 

Dividing  through  by  the  term  J~[  C,  and  adding  like  terms  yields: 

»=i 
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25  1  J_  1  >  + 

c  i  c2  c3  c4 


2  I  1  1  |  1  1  t  1  t  1 

1  1  c2c3  c2c4  c2c5  c3  c4  c3  c5  c4cs 


CiC3 

+ 

CjC4 

+ 

CiC5 

+ 

C3C4 

+ 

C3CS 

+ 

C4C5 

1 

1 

1 

A. 

1 

1 

4_ 

1 

C 1  c2 

“T* 

CiC4 

C1C5 

c2c4 

c2  c5 

c4c5 

1 

1 

1 

4. 

1 

4_ 

1 

4. 

1 

ctc2 

CiC3 

CiC5 

c2  c3 

c2  c5 

C3C5 

1 

1 

1 

1 

1 

1 

C\  c2 

"r 

CiC3 

CiC4 

c2  c3 

c2  c4 

C3C4 

|  C2  +  C3  +  C4  +  C  5  cx  +  c3  +  c4  +  c5 

“\2‘[  c2c3c4c5  J  +  Zj[  C,C3C4C5  j  + 

Ci  +  C2  +  C4  +  C5  c1+c2  +  c3  +  c5 

73  c1c2c4c5  J  +  *4  [  ClC2C3C5  J  + 

f.  ^  ' 

Cj  +  C2  +  C3  +  C4  z  j  z2 

5  C!C2C3C4  C2C3C4C5  CiC3C4C5 

z3  z4  z< 

- £ - + - - - + - - -  =  o 

CiC2C4C5  C*i  CTj  ^3  C5  C\C2C3C4 


(3.28) 


To  maintain  similarity  with  the  forms  of  the  quaternary  and  ternary  equations,  we 
can  separate  the  general  term  in  the  coefficient  for  a  in  Equation  (3.28)  into  four 
fractions: 


Cj  +  Ck  +  Cj  +  Cm  I  ^ 


Cj  Ck  Ci  Cm 


_ 1 _  _ 1 _  _ 1 _  1 

Zi[ckC,Cm  CjC'Cm  Cj  Ck  Cm  +  Cj  Ck  Q 

(3.29) 
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After  multiplying  each  of  the  five  fractions  by  zi  and  collecting  terms  with  common 
denominators,  the  following  form  appears: 


*.•  +  zj 

Cl  Cl  C„ 


(3.30) 


5 

We  can  invoke  £  zt  =1  to  construct  fractions  with  similar  subscripted  terms  in 

i=l 


both  numerator  and  denominator: 


\-zk-zl-zm 

Ck  Cl  Cm 


(3.31) 


This  yields  a  final  polynomial  expression  of  the  same  general  form  as  those  of  War¬ 


ren  (1991): 


a4  +  a3 


1  “  Z2  1  -  23  1  -  24  1  ~  25 

c.  +  c2  +  c3  +  c4  +  c5 


l-Zi~Z2  1  —  Zj  —  23  1  -  Z]  -  z4  1  —  2  1  ~  25  1  -  z2  -  Z3 

Ci  c2  cxcz  cx  c4  cx  c5  C2C3 


l-z2-zA  1  -22-z5 

c2  c4  C2C5 


1  -23-24  1  -  23  -  25 

C3C4  C3  C5 


1  -  z4  -  Z5 
c4  c5 


1  -2l  -Z2-Z3  1  —  2  1  ~  Z2  ~  24  1  -  21  -  22  -Z5  1  -  Zj  -  Z3  -  Z4 

CiC2C3  C,C2C4  C  XC  2C  $  CXC3C4 


1  -  zt  -  z3  -  25  1  -  -  z4  -  z5  1  -  z2  -  z3  -  z4  1  -  z2  -  z3  -  z5 
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This  yields  the  polynomial  expression  for  a  quinary  system: 


a4  +  a3 


\s  (1  -Zi) 

I  C- 

«=i 

+  a2 

y  3  0  ~z;> 

w  >=r+i  ci  Cj 


a 


4  5  6  (1  -  Zj  -  ZJ  -  Zj) 

Z*  2^  s*  r*  r* 

|_i=l  >=«'+!  k-i+2  Li 


5  zi 

I  V- 
i=1  n  C; 

j*i 


=  0  (3.33) 


3.3.2  Reduction  to  Quaternary  System 


Before  proceeding  to  develop  the  equations  for  six-  and  ^even-component  sys¬ 
tems,  we  must  ensure  that  the  quinary  equation  will  reduce  to  that  of  a  quaternary 
system  under  the  proper  conditions.  This  is  accomplished  by  setting  z5  equal  to 
zero  and  K5  equal  to  one  (Warren,  1991). 

When  z5  becomes  zero,  so  must  x5  and  y5.  This  would  seem  to  leave  K5 
undefined: 


lim  =  —  =  undefined  (3.34) 

Xj— »o  0 

ys-*o 

We  can  remove  this  difficulty  by  the  application  of  l’Hospital’s  Rule.  The  expres¬ 
sion  becomes: 


lim 

X,— >0 

y  3  *0 


<ty  5 
dy5 

& 5 
dx  5 


(3.35) 


Therefore,  C5  =  AT5  —  1  =  1  —  1  =  0. 
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5 

To  avoid  division  by  zero,  multiply  Equation  (3.32)  by  Jf  C,  : 

«= l 

ct4Cj  C2C3C4C5  + 

a3  [  (1  -  2  1)  C2  C3  C4  C5  +  (1  -  z2)  Cj  C3  C4  C5  +  (1  -  z3)  Cj  C2  C4  C5  + 

(1  -  z4)  C !  C2  C3  C5  +  (1  -  z5)  C  j  C2  C3  C4  ]  + 

a2  [  (1  -  *1  -  22)  C3  Q  C5  +  (1  -  zj  -  z3)  C2  C4  C5  +  (1  -  zj  -  Z4)  C2  C3  C5 

+  (1  -  2 1  — Z5)C2C3C4  +  (1  —  22  —  z3)  C 1  C4C5  +  (1  — z2  —  z4)C}  C3C5  + 
(1  ~  22  —  Z5)  Cj  C3  C4  +  (1  —  z3  —  z4)  Cj  C2  C5  +  (1  —  z3  —  Z5)  C\  C2  C4  + 

(1  -  z4-z5)C]  C2  C3]  + 

Ot  ^  Z  i  C  j  (C  2  +  C3  +  C4+C5)  +  Z2C2(Ci  +  C3  +  C4  +  C5)  + 

z 3  C 3  (C j  +  C2  +  C4  +  C5)  +  z4C4  (C 1  +  C2  +  C3  +  C5)  + 

2  5  C  5  (C 1  +  C2  +  C3  +  C4)j  +  z,  Ci  +  z2C2  +  z3C3  +  z4C4  +  Z5C5  =  0 

(3.36) 

Let  z5  and  C5  equal  zero: 
a3CiC2C3C4  + 

a2  [  (1  -  z,)  C2  C3  C4  +  (1  -  z2)  Cj  C3  C4  +  (1  -  z3)  Cj  C2  C4  + 

(1  -  Z4)  Cj  C2  C3  j  + 

a  [rjC,  (C2  +  C3  +  C4)  +  z2C2(Ci  +  C3  +  C4)  +  z3C3(Ci  +  C2  +  C4)  + 

z4C4(Cj  +  C2  +  C3)  +  Zj  Cj  +  z2C2  +  z3C3  +  z4C4  =  0  (3.37) 
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This  is  identical  to  Equation  (3.21)  in  WarTen’s  work,  which  is  the  expanded  form 
of  the  quaternary  equation.  We  can  now  safely  derive  the  expressions  for  six-  and 
seven-component  systems. 

3.3.3  Senary  System 


a5  +  a4 


6  (1  ~  Zj) 

h  ci 


+  or 


5  6  (1 

2—)  Ami  /->  f 

i=l  j=i+ 1  L'i 


<r 


4  5  6  (1  -  Z,  -  Z;  -  Zk) 

ZEE  r!7r  " 

i=l  >=i+ 1  *=i+2  W  W 


a 


3  4  5  6  (1  -  Zj  ~  Z;  ~  -  Z/) 

L2a  2a  2d  2a  c  c  c  c 

i= 1  7=j+1  *=i+2  /=»+3  W 


A  2, 


+  s  6 
n c, 

j*» 


=  0 


(3.38) 


3 .3.4  Septenary  System 


a6  +  a5 


7  (1  -  Zj) 

1-1 


+  or 


y  y  (1  ~  z«  ~  zjl 
k  Ci  Cj 


or 


5  6  7  (1  -  Zj  -  Z;  -  Zk) 

ZEE  ■  cVc  ~ 

i=l  y=i+i  *=i+2  W  '■'* 


or 


4  5  6  7  (1  -  Zt  -  Zj  ~  Zk-  Zj) 

E  E  E  E — c  CV>- 

i=l  y=i+l  k=i+2  1=4+3 


3  4  5  6  7  (1  -Zi-Zj-Zk-Zl~  Zm) 

2a  2d  2a  2a  2a  c  c  c  c  —  c 

;= M  fc=i+2  /= i+Z  m=i44  H  W 


+  i-r1- 

,=l  n  c, 

j** 


=  0 


(3.39) 
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3.4  Formulation  of  a  Generalized  Equation  for  the  Vapor  Fraction 


The  objective  function  can  be  recast  in  the  following  form: 

N  N 

I  zi  Q  II  +  a  Cy)  =  0,  where  C,,  C;  =  constant  (3.40) 

i  =1  j  *i 

Equation  (3.40)  is  in  the  form  of  the  generating  function  for  the  elementary  sym¬ 
metric  functions,  ar : 


17  (1  +  t  *,)  = 


i= l 


According  to  Macdonald  (1979),  a0(*i)  =  1  and  ar (xt )  =  0  for  all  r  >  n. 


(3.41) 


We  can  now  express  the  objective  function  in  terms  of  the  r-th  elementary  sym¬ 
metric  function  in  C,  : 

X  z,  C,  n  (1  +  a  Cj )  =  X  ci  X1  C, . CN)  =  0  (3.42) 

i=l  j*i  i=l  r= 0 


where  C,  indicates  the  exclusion  of  the  j'-th  term  from  the  operation. 
Since  ar  does  not  involve  i,  we  can  invert  the  order  of  the  summations: 


N- 1 


X  a' 


|  X  zi  Q  ar(Cl .  Q»—» 


=  o 


(3.43) 


A  working  definition  of  the  elementary  symmetric  function  ar  could  be  "tak¬ 
ing  permutations  of  the  elements  of  a  set  r  terms  at  a  time."  For  example, 


Qi(Ci,  Ctf)  =  (Cj  +  C  2+  +  Cn) 


(3.44) 
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The  condition  C(  is  equivalent  to  the  j*i  condition  imposed  on  the  summa¬ 
tion  terms  in  the  earlier  versions  of  the  a-polynomial  and  herein  lies  the  computa¬ 
tional  awkwardness.  We  want  to  find  an  expression  that  allows  the  summation  to 
proceed  over  all  N  components,  which  is  an  operation  readily  represented  by  a 
DO-loop  in  computer  programming. 

To  eliminate  Ci ,  we  must  expand  the  symmetric  function.  In  Chapter  4,  we 
will  tackle  this  problem  after  a  discussion  of  symmetric  functions. 


Chapter  4 


DEVELOPMENT  OF  THE  GENERALIZED  EQUATION 


In  this  chapter,  we  shall  present  a  brief  introduction  to  the  theory  of  symmetric 
functions  to  show  why  they  provide  such  a  powerful  tool  to  express  permutations. 
Then  we  will  show  the  reasoning  used  in  the  search  for  a  recursive  expression  for  a 
in  terms  of  N ,  Ci  and  z, .  Finally,  we  will  present  a  generalized  multicomponent 
equation  for  the  vapor  fraction,  a,  that  is  compact  and  readily  programmed  on  a 
computer. 

4.1  Introduction  to  Symmetric  Functions 

4.1.1  Notation  and  Definitions  of  Partitions 


Any  collection  of  v  non-negative  integers  (excluding  zero)  whose  sum  is  w  is 
called  a  v-partition  of  w.  The  individual  integers  are  referred  to  as  parts  of  the  par¬ 
tition  and  are  conventionally  written  in  descending  order  of  magnitude. 

David  et  al.  (1966)  state  that  if  there  are  X  distinct  parts,  say  pi,  p2,...,  p*. 
with  Pi>P2>Pi>  •  •  *  >Px  2:  1  and  if  p,  is  repeated  Jt,  times,  with  i  =  1,  2,...,  X, 

then  the  partition  is  written  (p  *'p  22  .  .  .  P\~).  The  weight,  w,  of  the  partition  is 
written  as 

x 

*  =  LA*. 

i=l 


(4.1) 
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and  the  number  of  pans,  or  length ,  is 

v  =  £  (4.2) 

i=i 

Macdonald  (1979)  refers  to  Tt,  as  the  multiplicity  of  /  in  the  partition.  For  example, 
the  partition  (42213)  has  weight  13,  6  pans  and  3  distinct  pans.  In  our  notation, 

P\  =  4  and  rtj  =  2;  p2  =  2  and  7^  =  1;  p3  =  1  and  rc3  =  3 
4.1.2  Symmetric  Functions 


A  symmetric  function  is  one  in  which  the  individual  pans  can  be  interchanged 
without  altering  the  value  of  the  function,  such  as 

=  xx  +jc2 +  *3+  •  •  •  +  *»  (4.3) 

i=i 

The  number,  n,  of  the  quantities  x  does  not  affect  the  relationships  between  the 
various  fomis  of  the  symmetric  functions,  but  does  appear  in  the  final  expressions. 
David  et  al.  (1966)  write 

X  x,  =  (1),  2  xi  ~  (r)  and  Z  xixj  ~  (")»  for  r  *  s  (4.4) 

i=l  i=l  i*j 

This  leads  directly  to  the  definitions  of  two  special  forms  of  symmetric  functions. 
MacMahon  (1920)  defines  the  unitary  or  a- functions  as 

a,  =  (  T  )  =  £  x.y.JCj,  ,  r  =  1,  2,  •  ••  (4.5) 

i  !<...<!, 

and  the  power  sums,  or  j-functions,  as 

=  (»■)  =  Z  xi  *  r  =  ’ 

i=i 


(4.6) 
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A  special  case  of  the  a-function  is  the  augmented  unitary  symmetric  function,  ur 
(David  et  al.,  1966): 

ur  =  [  lr  ]  =  r!(  V  )  =  r!ar  =  £  xil-xir  -  (4.7) 

summed  over  all  ordered  sets  ilf  .  .  .  ,  ir. 

4.1.3  Recursive  Expressions  for  Symmetric  Functions 

4.1.3.1  Interexpressibility  Tables 


Roe  (1918)  compiled  comprehensive  interexpressibility  tables  relating  the  vari¬ 
ous  classes  of  symmetric  functions  to  one  another.  These  consist  of  a  matrix  of 
coefficients  to  be  used  in  a  polynomial  which  might  yield,  for  example,  ur  =  /  (rr ). 
Of  interest  to  this  work  is  her  relationship  between  the  a-functions  (often  called  ele¬ 
mentary  functions)  and  the  s-functions;  it  is  presented  here  in  a  form  more  clearly 
expressed  by  David  et  al.  (1966): 


a,  = 


1  *  (_i  )<'+■*) 

r\  r  **  ** 


ni 

C  1  C  * 

*P  1 


m=l  7  *l!-*X*  pl 


■P\ 


(4.8) 


David  et  al.  (1966)  also  used  this  equation  to  construct  interexpressibility 
tables  describing  polynomials  in  power-sum  series  (s)  for  a-functions  up  to  and 
including  weights  of  12.  For  instance,  a  unitary  symmetric  function  of  weight  3 
would  be  represented  by  the  following  polynomial  from  their  Table  1.5.3: 

°3  =  3r[(1)3-3{2)(1)  +  2(3)] 


(4.9) 


39 


which,  in  terms  of  s-functions,  is: 

fl3=— [^j3  -  3s2s  1  +  2j3  j 

and,  when  written  as  power  sums,  becomes: 


(4.10) 


a3  =  ¥ 


' 

•  -N 

n 

Z 

i=i 

"  .2  " 


-  3X  xi  Z  xi  +  2Z  xi 

i= 1  i= 1  1=1 


(4.11) 


However,  neither  Equation  (4.8)  nor  Equation  (4.11)  is  conducive  to  solution  by 
computer  without  a  tremendous  table  look-up  effort. 


4.1.3.2  Determinant  Form 


Fortunately,  David  et  al.  (1966)  present  another  relationship  between  ar  and 
sr  in  determinant  form: 


*i 

1 

0 

0  • 

•  0 

s2 

*1 

2 

0  • 

•  0 

53 

*2 

3  • 

•  0 

det  M 

S* 

*3 

52 

*i  • 

•  0 

r ! 

• 

• 

• 

. 

•  r-1 

*r 

*r-l 

S,  -2 

*r-3  • 

•  *i 

(4.12) 


This  provides  a  practical  method  of  calculating  ar  that  is  also  readily  programm¬ 
able  on  a  computer. 
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4.2  Search  for  a  Recursive  Expression  for  the  Vapor  Fraction 

Armed  with  a  working  knowledge  of  symmetric  functions,  we  can  manipulate 
the  a  equation  developed  in  Chapter  3: 


N- 1 

X  «r 


r= 0 


N 


X  zi  Q  ar(C  1,...,  Cj,...,  Cyv) 


<=1 


=  0 


(3.43) 


to  eliminate  the  exclusion  term  C,  and  expand  the  symmetric  functions  into  a  more 
recognizable  form.  We  will  examine  the  results  for  several  values  of  r  and  use 
them  to  write  a  general  expression  for  a  as  a  function  of  N. 


4.2.1  Case  I:  r  =  N-2 


Equation  (3.43)  yields  the  following  coefficient  for  a: 


,<N-2) 


N 

X  zi  C'l  aN-l(C  lvi  Cj  »•••»  CN) 

i=l 


(4.13) 


We  can  expand  the  symmetric  function  aN_2  as  shown  in  Equation  (4.14).  Since 
the  exclusion  of  C,  from  the  product  on  the  RHS  gives  (ZV-1)  terms,  we  must  sub¬ 
tract  -pr  from  the  sum  to  yield  (N-2): 


aN-2 (C !»...,  Cj,...,  C^)  —  (Cr  Cr-CN) 


■4-  +  -2~+ 

C 1  c2 


(4.14) 


To  eliminate  C; ,  we  can  write  the  product  on  the  RHS  of  Equation  (4.14)  as 


(Cx  -Ci -CN)  = 


(Ci"Cy) 

C. 


(4.15) 


This  maneuver  will  allow  the  summation  to  proceed  over  all  N  components. 
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After  substituting  Equation  (4.15)  into  Equation  (4.14),  we  have: 


(CfCjv) 

JL  + 

_L+  .. 

i 

1 

C, 

c, 

c2 

CN 

"  c( 

(4.16) 


Substituting  Equation  (4.16)  into  Equation  (4.13),  cancelling  C, ,  multiplying  by  z, , 
and  then  summing  over  i  gives: 


a' 


(N-2)  J 


(C,~CW) 


rN  ' 

E  z« 


Cl  C2 


>  M  1 

N  z- 

-E~ 

► 

*- 1  Ci 

m 

J 

(4.17) 


A/  A/ 

We  recall  that  £  Z;  =  1  and  recognize  that  (C^-C/y)  =  JI  C*.  Noting  the  pres¬ 
et  *=i 

ence  of  an  elementary  symmetric  function  in 


+  ••*+_ 
Ci  CN 


,  we  can  write 


Equation  (4.17)  as: 
a(N-2). 


N 

nc* 

k-l 


1 

cN 


N  H 

-if 

i=i  Li 


(4.18) 


4.2.2  Case  H:  r  =  N-3 


Equation  (3.43)  now  becomes: 


a(*-3)J 


N 

^  z4-  C,-  Of^—2  (Ci,...,  C,-,...,  C„)|> 

i=i 


(4.19) 


We  can  expand  the  symmetric  function  aN_3  as  shown  in  Equation  (4.20).  We 

eliminate  C,  in  the  same  manner  as  in  Equation  (4.15)  and  remove  in  a  similar 

w 

1 


fashion.  But  this  also  deletes  the  term 


c,2 


which  is  necessary  to  cancel  the 
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corresponding  term  in  the  product.  Therefore,  we  must  compensate  by  adding  — : 

C,2 


on-3(C C,,...,  CN)  - 


(CrcN) 


Ci 


1  +^-+ 


C  i C  2  C  jC  2 


Cn-\Cn 


'AT 


Ci2 


(4.20) 


Substituting  Equation  (4.20)  into  Equation  (4.19)  and  making  consolidations  similar 
to  the  previous  development  yields: 


a 


CIV-3) 


N 

nck 

i 


a2 


1  1 

Cl’"’  Cyv 


-a  i 


1  1 
Ci . CN 


*  Z:  N  Z: 

y  —  +  y  — 

^  r  "  /-2 

i=l  i=i  c, 


(4.21) 


4.2.3  Case  HI:  r  =  0 


We  have  saved  consideration  of  this  case  for  last  because  the  properties  of  a0 
are  not  readily  apparent.  It  would  seem  reasonable  to  interpret  a0(C1,...,  C, ,...,  CN) 
as  meaning  "taking  permutations  of  the  elements  of  a  set  zero  terms  at  a  time." 
However,  when  r  =  0,  ar  -»1  and  we  know  from  previous  developments  that  our 
a-polynomial  does  have  a  constant  term.  Therefore,  a0(Ci,...,  C,,...,  CN)  must 
equal  one ,  after  Macdonald  (1979).  So,  for  r  =  0,  Equation  (3.44)  becomes: 

r  \ 

N 

-  Z  2ici  ' 

i= 1 


(4.24) 
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We  show  Equation  (4.18)  and  Equation  (4.21)  again  to  look  for  patterns  that  may 
assist  us  in  writing  the  expression  for  the  (N-p)-\h  term: 


After  substituting  for  r ,  Equation  (3.43)  becomes: 


£  aF-fi)  £  2i  Ci  %^(C„...,  C, CN)\  =  0 
P=  1  i=  1  J 


(4.23) 


By  continuing  the  expansion  of  this  equation  in  the  same  fashion  as  in  the  first  two 
cases,  we  note  a  descending  order  of  the  symmetric  function  and  an  ascending 
exponent  of  C,  with  each  additional  term.  This  leads  to  a  general  expression: 

£  o^>  -  £  Zj  (Cr-c^L^Cf1 . c#1)  - 

p= i  i=i  L 

Crxap_2(C f1 ,...,  C/71)  +  C,~2a/,_3(Cf1 ,...,  C^1)  -  C~^ap_^(C  f1 ,...,  C^-1) 


+  •••  ±C1-^-2)a1(Cr1,...,C/71)  tCf^  1=0 


(4.24) 
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N 

Multiply  by  zt  ,  sum  i  from  1  to  N  and  recall  that  (C  fQv)  -nct: 

k- 1 


N 

£  a{N~p) 
p= i 


'  n 

t=i 

V 


N  ii 

—  Op-2  S  +  ap- 3 

i=l  W 


I 


i=l 


^  2 
'i 


+ 


£1 


=  o 


(4.25) 


N 

Since  n  Ck  does  not  involve  p,  we  can  move  this  term  outside  the  summa- 
*=i 

tion  sign  and  then  divide  it  out  as  a  factor  common  to  all  powers  of  a.  By  examin¬ 
ing  the  relationship  between  p,  the  subscripts  of  a  and  the  superscripts  of  C,  ,  we 
can  collapse  Equation  (4.25)  into  a  more  compact  form: 


where  ap_j  =  ap_}(  Cf1 ,...,  CNl  ) 


a  0=  1 

Ci  =  (Ki  )t,  p  ~  1 


(4.26) 

(4.27) 

(4.28) 

(4.29) 


Using  the  determinant  expression  for  the  elementary  symmetric  functions  that  was 
presented  in  4.1.3.2,  Equation  (4.27)  becomes: 


,  n  r-\\-  det  M 

i  . c*  >  (p ~j )! 

The  matrix  M  has  dimensions  (p-j)  x  (p-j)  and  elements  given  by: 


K/l  =i 


Sk-i+ i  if  /  <  k 
k  if  /  =  ifc+l 
0  if  /  >  ifc+l 


N 

The  s  elements  are  given  by  s*  =  £ 

i= l 


f  ^ 


,  X*  1,  2 . (p-j) 


(4.30) 


(4.31) 


C, 


(4.32) 
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VALIDATION  OF  THE  GENERALIZED  EQUATION 


The  first  test  of  validity  for  Equation  (4.26)  requires  that  it  be  equivalent  to  the 
form  of  the  objective  function  presented  in  Equation  (3.40).  Second,  it  must  gen¬ 
erate  the  same  coefficients  for  the  a  polynomial  that  were  produced  through  the 
expansion  of  the  objective  function  in  Equation  (3.25)  through  Equation  (3.32). 
Third,  the  equation  must  predict  the  proper  vapor  fraction  for  a  fluid  undergoing  an 
isothermal,  isobaric  flash  process. 

The  first  test  is  supplied  by  a  mathematical  proof  in  Appendix  A.  The  second 
test  can  be  accomplished  by  comparing  the  coefficients  produced  by  Equation  (4.26) 
with  those  of  Equation  (3.33).  Since  this  equation  has  already  been  shown  to 
reduce  to  that  for  a  quaternary  system  under  the  proper  constraints  on  z5  and  K5, 
then,  by  induction,  we  can  state  that  the  polynomial  produced  by  Equation  (4.26) 
will  do  the  same  and  therefore  should  be  valid  for  any  number  of  components. 

The  third  test  will  be  satisfied  by  comparing  the  equilibrium  ratios  generated 
by  Equation  (4.26)  with  experimental  values  determined  for  several  multicomponent 
hydrocarbon  fluids. 

5.1  The  Generalized  a  Equation  for  a  Quinary  System 


For  a  five-component  system.  Equation  (4.26)  becomes: 


(5.1) 
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which  will  yield  a  quartic  polynomial  in  a: 


p.4a4  +  jjL3a3  +  )j.2a2  +  jija  +  }Xq  =  0 


5.1.1  Coefficient  (p  =  1) 


U4  =  fl0(Q_1)  [  z !  4-  z2  +  2 2  +  z4  +  Z5  ]  (5.3) 

We  have  already  said  that  a0(C,-1)  is  defined  as  one  and  the  sum  of  the  mole  frac¬ 
tions  also  equals  one,  so  Equation  (5.3)  yields: 


M-4  =  1 


i=l  «=1  W 


*■>-  i+i+i+i+i  (1)- 


Z\  z2  z3  *4  z5 

<»  c7+c7+c7+c7  +  ^- 


1-Zj  1  —  Z2  1  —  Z3  1  —  Z4  1  —  Z5 

=  c[  +  — +  c3  +  c4  +  c5 


(5.7) 
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5.1.3  Coefficient  lb  (p  =  3) 


H2  =  a2(Cf'>  £  *i  -  ai(Cr')  £  7T-  +  «o<Cf'>  £  77 

i=l  i=l  i=l  W 


(5.8) 


M-2  = 


1  +^-  +  ^-  +  ^  +  ^L-  +  ^  +  ^+  1 


C,c2  C,C3  C1C4  C,C5  C2C3  C2C4  C2C5  C3C4 


1  1 
+  —  + 


C3C5  C4C5 


(D- 


11111 
Cj  c2  c3  c4  c5 


il  il  il 

Ci  c2  c3 


+  (1) 


+iL  +  i2.  +  ii.  +  iL 

‘  c42  c52 


c,2 


c } 


cl 


(5.9) 


^2=  CiC2  +  C1C3  +  C,C4  +  C1C5  +  C2C3  +  C2C4  +  C2C5  +  C3C4  + 


1 _  1  _  Jj _ z_2 _ Z3 _ ZJ_  _  Z1  +  Z2  _  z\  +  *3 

C3C5  +  C4C5  c2  cl  cl  Cl '  c|  C1C2  CiC3 


Z  i  +  Z4  Zi  +  25  Z2  +  z4  Z2  4-  Z5  Z3  +  Z4  Z3  +  Z5 

C1C4  C1C5  C2C3  C2C4  C2C5  C3C4  C3C3 


n2  = 


z4  z5  Z1  z2  z3  z4  Z5 

C4C5  +  C,2  C22  C32  C42  C52 

1  -  Zi  -  Z2  1  ~  Zj  ~  Z3  1  -  Zj  -  z4  1  -  Zi  -  Z5 

CiC2  CiC3  C1C4  C1C5 


(5.10) 


1  -  z2  -  Z3 

C2C3  + 


1  -  z2  -  z4  1  -  z2  -  z5  1  -  Z3  -  z4  1  -  Z3  -  Z5  1  -  Z4  -  z5 
C  2C  4  C  2C  5  C  3C  4  C  3C  5  C  4C  5 


(5.11) 
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5.1.4  Coefficient  u.n  (p  =  4) 


Ml  =  t,  -  i2<C-')i  jr  +  0,(Cf>)£  -  ^(Cf1)  £  (5.12) 

1=1  1=1  1=1  Cj  i=i  c,- 


5  z, 


Mi  = 


1  +^1  +^_  +  ^J^+  >-  + 


C1C2C3  ^  1^-  2^4  CiC2C5  C 1C3C4  CJC3C5 


C1C4C5  C2C3C4  C2C3C5  C2C4C5  C3C4C5  J 

11111111 
C1C2  CjC3  C,C4  C2C3  C2C4  C2C5  C3C4 


1  ,  1  ' 

zi  z2  z3  z4  z5 

C3C5  C4C5 

C!  c2  c3  C4  C5 

_L  + JL  +  _L  + J_  + J_ 

2L  +  iL  +  iL  +  ii_  +  iL 

Ci  C2  C3  C4  C5 

2  2  /-» 2  2  /"*  2 
t  j  Cj  1 3  C  4  1 5 

(1) 


Zl  +lL  +  il  +  ii.+  *5 


c? 


c,3  c,3  cl  c<3 


(5.13) 


It  is  evident  that  the 


C?C, 


terms  in  the  second  part  of  Equation  (5.13)  will  cancel 


those  in  the  third  part,  while  the  — -  terms  in  the  third  part  will  negate  the  entire 

Q 

fourth  part  of  the  equation.  The  first  and  second  parts  yield: 


1  -  *t  -  *2-*3  1  ~*t  -  *2  -  *4  1  -zx-z2~z5  1  -  zi  -z3  -  z4 

CJC2C3  ^*1^2^ 4  C jC  3C  4 
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+  \  -  ZX-  Z  3-z5+  1  -  ZX  -  Z4  -  Z  5  +  1  ~  z2  ~  z3  ~  z4  +  1  ~  z2  ~  z3  ~  ZS  + 


C1C3C5 


C 1C4C5 


C  2C3C  4 


C  2C3C5 


1  -  z2  -  Z4  -  Z_5_  +  1  -  z3  “  z4  “  z5 


C  2C  4C  5 


C3C4C5 


(5.14) 


5.1.5  Coefficient  Uo  (P  =  5) 


Mo  =  fl«(C,-1)  i  »i  -  a3(C,-')  £7-  +  arfCf ')  £  77  - 
(=1  i=i  w  ;=i  w 


«i(C,-')  £^t  +  oo(C,'‘)  £  77 

i=l  i=l  W 


(5.15) 


The  analogous  cancellations  of  the  higher-order  -7-  terms  will  occur,  leaving  a 

w 


sum  of  five  terms  having  the  form 

1  ~  z,  ~  zj  ~  zk  ~  zi 

C,  Cj  Ck.  Q 


(5.16) 


Since  the  mole  fractions  must  sum  to  one,  we  can  replace  the  numerator  of  Equa¬ 
tion  (5.16)  with  the  mole  fraction  of  the  remaining  component,  zm  ,  to  yield: 


Ho  = 


C  2C  3C  4C  5 


C1C3C4C5 


C  \C  2C  4C  5 


z4 


C1C2C3C5  C1C2C3C4 


(5.17) 


A  term-by-term  comparison  with  Equation  (3.32)  shows  that  the  generalized  a  poly¬ 
nomial  [T  Ration  (4.26)]  produces  identical  results. 
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5 2  Reproduction  of  Experimental  Vapor-Liquid  Equilibrium  Data 

5.2.1  Flash  Calculation  Package 

The  flash  calculation  package  used  in  this  work  incorporates  the  /f-value  equa¬ 
tion  of  Wilson  (1969)  and  the  modified  PREOS  proposed  by  Stryjek  and  Vera 
(1986a).  An  attempt  was  made  to  use  the  Af-value  prediction  of  Varotsis  (1989) 
but,  as  noted  in  Chapter  2,  it  was  developed  to  characterize  a  broad- spectrum 
petroleum  reservoir  condensate  or  crude  oil.  It  experiences  difficulty  handling  an 
arbitrary  hydrocarbon  mixture,  such  as  the  artificial  systems  for  which  equilibrium 
data  is  available  to  validate  this  work. 

The  volume  correction  of  Mathias  et  al.  (1989)  and  the  complementary  calcu¬ 
lation  of  Schick  and  Prausnitz  (1968)  for  mixture  pseudo-critical  volume  are  incor¬ 
porated  into  the  PRSV  EOS  but  since  it  is  only  required  to  generate  compressibility 
factors  and  fugacities,  the  modifications  have  no  noticeable  effect  on  the  computa¬ 
tions.  The  PRSV  EOS  shows  marked  improvement  over  the  PREOS  when  used  to 
duplicate  bubblepoint  and  dewpoint  studies  performed  by  Warren  (1991). 

The  binary  interaction  coefficients  used  in  the  PRSV  EOS  are  taken  from 
Knapp  et  al.  (1982)  and  Walas  (1985).  Physical  property  data  and  equation  param¬ 
eters  for  the  chemical  components  are  extracted  from  Stryjek  and  Vera  (1986b,c), 
Kumar  (1987>  and  Proust  and  Vera  (1989). 

The  computation  of  the  determinant  used  to  generate  the  elementary  symmetric 
functions  is  accomplished  with  a  modified  Gaussian  elimination  routine.  The  first 
elementary  symmetric  function,  a\,  is  defined  by  a  [lxl]  matrix,  whose  deter¬ 
minant  constitutes  the  element  itself.  By  definition,  a0  is  set  equal  to  one. 
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The  polynomial  is  evaluated  at  the  bubblepoint  line  (a  =  0)  and  an  interval- 
halving  technique  is  used  to  inarch  across  the  two-phase  region  until  the  value  of 
the  polynomial  changes  sign,  indicating  the  vicinity  of  the  root.  Then  a  Newton- 
Raphson  iterative  search  is  conducted  to  converge  to  the  exact  value  of  a. 

5.2.2  Binary  System 

The  fugacity-based  flash  algorithm  is  used  to  replicate  the  equilibrium  ratios 
determined  by  Bloomer  et  al.  (1953)  for  a  methane-ethane  system  at  a  temperature 
of  -60  °F  over  a  pressure  range  of  100-900  psia.  Comparisons  of  calculated  and 
empirical  values  of  and  ^c2h6  appear  in  Figure  5.1  and  5.2,  respectively.  The 
results  lie  within  the  margin  of  error  attributable  to  the  PRSV  EOS. 

5.2.3  Septenary  System 

Standing  (1977)  provides  a  sample  flash  calculation  for  a  seven-component 
hydrocarbon  system  reported  by  Dodson  and  Standing  (1941),  complete  with  values 
for  experimental  Ki  and  the  vapor  fraction.  This  sort  of  data  allows  the  calculation 
of  a  solely  on  the  basis  of  computing  the  coefficients  of  the  a-polynomial  and 
determining  the  applicable  root,  with  no  recourse  to  the  equation  of  state.  Once  the 
interval-halving  search  provides  an  initial  estimate  of  the  root,  the  Newton-Raphson 
technique  converges  in  three  iterations  to  a  value  of  a  identical  to  that  calculated  by 
Standing. 


Equilibri 
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5.2.4  Predicting  Roots  with  the  Fourier-Budan  Theorem 

A  useful  theorem  for  predicting  the  number  of  roots  of  a  polynomial  that  can 
occur  on  a  particular  interval  is  that  of  Fourier  and  Budan  (Barbeau,  1989).  Sup¬ 
pose  p(t)  is  a  polynomial  over  the  field  of  real  numbers,  R,  and  that  u  and  v  are  real 
numbers  with  u  <  v  and  p(u)p(v )  *  0 .  The  number  of  zeros  between  u  and  v 
cannot  be  greater  than  A  -  B ,  where  A  is  the  number  of  changes  of  sign  in  the 
sequence  {  p(u),  p'(u),  p"(u),  ....  p^n\u)  }  and  B  is  the  number  of  changes  of 
sign  in  the  sequence  {  p  (v),  p'(v),  p"(v p(n)(v)  }•  If  this  number  differs 
from  A  -  B ,  it  must  do  so  by  an  even  amount. 

An  interesting  aspect  of  the  polynomial  expression  for  the  vapor  fraction  is  its 
capability  to  mathematically  confirm  the  existence  of  a  unique  value  within  the 
two-phase  region  for  a  particular  set  of  feed  conditions.  This  is  equivalent  to  stat¬ 
ing  that  the  polynomial  has  only  one  zero  on  the  interval  0  <  a  <  1.  From  the  phy¬ 
sics  of  the  problem,  we  know  this  to  be  true  but,  by  the  use  of  the  Fourier-Budan 
theorem,  we  can  also  prove  it  rigorously. 

Let  us  test  this  theorem  on  the  septenary  system  of  Standing  (1977)  utilized  in 
5.2.3;  this  is  represented  by  a  sixth-order  polynomial: 

ma6  +  it5a5  +  p4a4  +  p3a3  +  \i2a2  +  pja  +  ^  =  0  (5. 18) 

where 


Ho  = 

-9.58519 

H3  = 

87.24949 

Hi  = 

65.90501 

H4  = 

-21.71701 

H2  = 

-120.72959 

Hs  = 

-1.76522 

m  =  1.00000 


56 


We  can  differentiate  Equation  (5.18)  six  times  and  form  the  derivative  sequences 
for  u  =  0  and  v  =  1.  The  sign  changes  are  summarized  in  Table  5.1. 


Table  5.1  -  Derivative  Series  of  Fourier-Budan  Theorem: 
7-Component  Hydrocarbon  System  (Standing,  1977) 
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Since  A  -  B  =1,  there  exists  only  a  single  root  of  the  polynomial  on  the 
interval  [0,1];  an  examination  of  the  graph  of  the  function  (Figure  5.3)  confirms  this 
fact.  Therefore,  we  can  use  the  interval-halving  and  Newton-Raphson  solution  pro¬ 
cedures  outlined  at  the  beginning  of  this  chapter  with  confidence  that  they  will 
obtain  a  unique,  realistic  value  of  the  vapor  fraction. 

5.2.5  Decenary  System 

Gregory  et  al.  (1971)  performed  equilibrium  measurements  on  a  lean  natural 
gas  at  cryogenic  conditions.  It  is  reported  as  a  ten-component  system  with  the  feed 
composition  shown  in  Table  5.2.  This  is  a  very  "sparse"  ten-component  gas,  with 
six  components  present  in  extremely  dilute  concentrations.  The  /((-values  for  the 
last  four  constituents  were  zero  for  eleven  of  the  sixteen  operating  conditions  tested 
in  this  work,  denoted  in  Table  5.3  by  the  run  number  assigned  by  the  investigators. 
The  remaining  twelve  sets  of  published  data  duplicate  conditions  in  one  of  the 
tested  runs  or  are  incomplete  due  to  apparatus  failure.  The  use  of  the  Fourier- 
Budan  theorem  provides  warning  that  perhaps  this  gas  would  be  better  represented 
by  an  equivalent  "lumped"  system. 

Recall  that  the  number  of  roots  predicted  by  the  Fourier-Budan  theorem  is  the 
maximum  possible  and  may  differ  from  the  true  value  by  only  an  even  integer. 
This  is  demonstrated  in  Table  5.4,  where  both  the  predicted  and  actual  number  of 
roots  for  each  run  are  tabulated.  The  Newton-Raphson  technique  converges  to  the 
experimental  value  for  ten  of  the  sixteen  runs.  Three  other  data  points  follow  the 
proper  trend,  while  no  root  is  found  on  the  interval  [0,1]  for  three  other  conditions 
(Figure  5.4). 


Table  5.2  •  Feed  Composition: 
10-Component  Natural  Gas 
(Gregory  et  ai.,  1971) 

Component 

Component 

h 

Nitrogen 

0.00600 

n-Butane 

0.00070 

Methane 

0.95970 

/-Pentane 

0.00030 

Ethane 

0.03000 

n-Pentane 

0.00010 

Propane 

0.00390 

3-Methylpentane 

0.00025 

i-Butane 

0.00070 

2-Methylhexane 

0.00015 

Table  5.3  -  Experimental  Flash  Conditions: 
10-Component  Natural  Gas 
(Gregory  et  alM  1971) 


Run 

Pressure 

(psia) 

Temperature 

(°F) 

Run 

Pressure 

(psia) 

Temperature 

(°F) 

1 

300.0 

-156.3 

14 

100.0 

-200.0 

3 

100.0 

-206.0 

15 

500.6 

-127.0 

D 

700.0 

-103.0 

18 

23.0 

-252.0 

B 

500.0 

-125.0 

20 

497.0 

-129.0 

8 

498.5 

-120.0 

21 

23.5 

-251.5 

9 

695.0 

-105.0 

25 

700.0 

-107.0 

10 

100.0 

-203.3 

26 

298.0 

-157.5 

12 

100.0 

-195.0 

28 

500.0 

-130.0 

Table  5.4  -  Results  of  the  a-Polynomial  and  Fourier-Budan  Theorem: 
10-Component  Natural  Gas 
(Gregory  et  at.,  1971) 


Run  Root  Limit 
(Actual) 


Newton- 

Raphson 

Iterations 


Roots  on  [0,1] 


Initial  Guess 

Calculated 

0.605 

0.603 

0.825 

0.822 

***** 

0.775 

0.772 

0.905 

0.904 

***** 

***** 

0.695 

0.692 

0.915 

0.912 

0.895 

0.891 

0.965 

0.966 

0.835 

0.830 

0.735 

0.737 

0.845 

0.843 

0.045 

0.044 

0.415 

0.415 

0.585 

0.587 

0.775 

0.773 

0.015 

0.011 

0.645 

0.642 

0.385 

0.380 

0.435 

0.434 

0.485 

0.480 

0.155 


0.911 


0.761 


0.908 


0.795 


0.830 
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An  examination  of  graphs  of  the  polynomial’s  behavior  over  a  range  of  a  for 
Runs  3  and  4  (Figures  5.5  and  5.6)  confirms  the  algorithm’s  prediction  that  no  roots 
exist  within  the  phase  envelope.  The  case  of  Run  9  is  not  so  obvious.  Its  graph 
(Figure  5.7)  shows  that  the  function  exists  entirely  above  the  abscissa;  hence  no 
root  is  possible.  However,  if  the  resolution  of  the  graph  is  increased  to  examine  the 
region  very  near  the  axis,  two  local  minima  are  revealed  (Figure  5.8).  One  of  these 
corresponds  to  the  experimental  value  of  a  determined  for  this  run.  The  polyno¬ 
mial  is  attempting  to  represent  the  system’s  behavior  but  is  not  completely  success¬ 
ful  because  the  low  concentration  of  certain  components  effectively  prevents  the  gas 
from  acting  like  a  decenary  system. 

It  is  instructive  to  compare  the  form  of  the  a-polynomial  with  that  of  the 
Rachford-Rice  objective  function  which  is  superimposed  on  Figure  5.7.  The  same 
high-resolution  scan  of  the  graph  of  the  latter  equation  depicts  no  equivalent  max¬ 
ima  which  might  identify  the  vapor  fraction  in  the  manner  of  the  polynomial. 

5.2.6  Lumping  a  Decenary  System  into  a  Quaternary  System 

The  a-polynomial  successfully  converges  to  the  proper  answer  for  a  majority 
of  the  runs;  however,  it  also  yields  multiple  roots  where  the  physics  of  the  problem 
allows  only  one.  This  suggests  that  the  system  is  not  being  properly  modeled.  The 
categorization  of  the  fluid  as  a  ten-component  natural  gas  is  overly  generous  in  light 
of  the  fact  that  six  of  its  chemical  constituents  are  present  in  mole  fractions  meas¬ 
ured  in  the  ten-thousandths.  It  was  decided  to  represent  this  sparse  gas  as  a  four- 
component  lumped  system,  consisting  of  methane,  ethane,  nitrogen  and  propane. 

The  mole  fractions  of  this  new  fluid  are  normalized  and  the  resulting  cubic 
polynomial  in  a  is  solved.  The  Fourier-Budan  theorem  predicts  a  maximum  of  one 
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root  in  the  two-phase  region,  to  which  all  sixteen  runs  converge.  The  numbers 
tabulated  in  Table  5.5  and  displayed  graphically  in  Figure  5.9  attest  to  the  validity 
of  this  lumping  scheme.  An  attempt  was  made  to  eliminate  the  next  leanest 
component-propane-ffom  the  mixture  and  model  the  system  as  a  ternary,  but  this 
resulted  in  spurious  roots  for  all  data  runs  and  was  hence  rejected  as  unrealistic. 


Table  5.5  -  Results  of  the  a-Polynomial  and  Fourier-Budan  Theorem: 
"Lumped"  4-Component  Natural  Gas 
(Gregory  et  ah,  1971) 

Run 

Root  Limit 

Newton- 

Raphson 

Iterations 

Roots  on  [0,1] 

(Actual) 

Initial  Guess 

Calculated 

Experimental 

1 

1 

3 

0.635 

0.631 

0.603 

3 

1 

3 

0.035 

0.038 

0.155 

a 

1 

3 

0.995 

0.998 

C  911 

D 

1 

3 

0.815 

0.819 

0.761 

8 

1 

3 

0.935 

0.932 

0.908 

9 

1 

3 

0.925 

0.920 

0.795 

10 

1 

3 

0.715 

0.711 

0.687 

12 

1 

3 

0.895 

0.899 

0.890 

14 

1 

3 

0.845 

0.842 

0.830 

15 

1 

3 

0.765 

0.770 

0.747 

18 

1 

3 

0.035 

0.037 

0.109 

20 

1 

3 

0.635 

0.631 

0.591 

21 

1 

3 

0.015 

0.012 

0.078 

25 

1 

3 

0.775 

0.779 

0.548 

26 

1 

3 

0.475 

0.473 

28 

1 

3 

0.535 

0.536 

0.486 

Chapter  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Conclusions 

1.  The  Rachford-Rice  objective  function  can  be  represented  as  a  polynomial  in  a, 
the  system  vapor  fraction.  Its  coefficients  involve  elementary  symmetric  func¬ 
tions,  which  can  be  expressed  in  terms  of  a  determinant  whose  elements  are 
functions  of  equilibrium  ratios  and  feed  composition. 

2.  The  polynomial  has  been  shown  to  be  well-behaved  within  the  two-phase 
vapor-liquid  region  if  the  system  is  properly  defined  in  terms  of  the  number  of 
its  components.  The  vapor  fraction  root  on  the  interval  [0,1]  can  be  quickly 
determined  using  an  ordinary  interval-halving  technique  to  provide  an  initial 
estimate  to  the  Newton-Raphson  iterative  method. 

3.  The  regular  behavior  of  the  polynomial  lends  itself  to  use  as  a  descriptive  tool 
for  the  conditions  of  the  system  within  the  phase  envelope.  The  Rachford- 
Rice  objective  function  is  not  capable  of  this  task  as  evidenced  by  Figure  5.7; 
its  unpredictable,  singular  nature  offers  no  clue  to  the  reason  a  root  was  not 
found  on  the  interval  [0,1]  for  this  case.  As  discussed  earlier,  a  close  exami¬ 
nation  of  the  curve  of  the  polynomial  revealed  a  local  minimum  at  the  experi¬ 
mental  value  of  a.  This  became  a  realistic  root  (a  <  1)  once  the  system  was 
lumped  into  four  components. 
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4.  The  theory  of  polynomials  is  well-developed  and  their  behavior  and  zeros  can 
be  predicted  with  good  confidence.  By  the  use  of  the  Fourier-Budan  theorem, 
it  can  be  shown  mathematically  that  only  one  real  root  for  the  a-polynomial 
can  exist  on  the  interval  [0,1]  for  a  well-defined  system.  This  eliminates  the 
need  to  solve  for  all  the  roots  of  an  jV-th  order  polynomial. 

5.  The  Fourier-Budan  theorem  can  be  used  as  a  tool  for  investigating  various 
lumping  schemes  whereby  multicomponent  fluids  are  condensed  to  equivalent 
systems  composed  of  fewer  components.  The  phase  behavior  of  sparse  fluids 
having  dilute  concentrations  of  several  constituents  does  not  seem  to  be  well- 
described  by  the  polynomial  of  degree  appropriate  to  the  number  of  com¬ 
ponents.  In  this  case,  the  polynomial  yields  no  roots  or  at  least  two  roots 
inside  the  phase  envelope  for  certain  temperature  and  pressure  conditions.  It 
appears  that  a  lumping  scheme  can  be  tuned  by  generating  pseudocomponents 
to  give  successive  polynomials  of  lower  degree  until  only  one  root  is  deter¬ 
mined  on  die  interval  0  <  a  <  1 . 

6.2  Recommendations 

1.  Further  study  should  focus  on  coupling  the  polynomial  algorithm  to  an  equa¬ 
tion  of  state  and  extending  this  work  to  flash  calculations  involving  more  than 
two  phases. 

2.  Timing  studies  could  be  conducted  to  determine  the  exact  savings  in  CPU  time 
realized  by  the  use  of  the  polynomial  instead  of  the  Rachford-Rice  objective 
function. 
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3.  Peng  et  al.  (1975)  estimate  that  75%  of  the  total  computing  time  in  composi¬ 
tional  reservoir  simulation  may  be  related  to  the  phase-behavior  part  of  the 
program.  The  savings  in  computational  workload  realized  by  the  use  of  the 
generalized  equation  developed  in  this  work  might  be  applied  to  the  employ¬ 
ment  of  an  EOS  better  able  to  predict  fluid  thermodynamic  properties.  The 
highly  nonlinear  nature  of  the  equations  proposed  by  Benedict,  Webb  and 
Rubin  (1940,  1942,  1951)  or  Lee  and  Kessler  (1975)  require  iterative  solutions 
but  they  yield  much  more  accurate  representations  of  fluid  behavior,  espe¬ 
cially  of  nonhydrocarbon  systems. 

4.  Since  the  coefficients  of  the  generalized  polynomial  depend  only  on  the  feed 
composition  and  equilibrium  ratios,  research  should  continue  to  develop  highly 
accurate  K-value  prediction  methods  (e.g.,  on  the  basis  of  convergence  pres¬ 
sure).  If  this  can  be  done  with  sufficient  accuracy,  the  fugacity-convergence 
approach  and  its  inherent  dependence  on  an  equation  of  state  can  be  sup¬ 
planted  for  flash  calculations  where  nothing  more  than  the  phase  split  and 
compositions  are  required.  The  polynomial  algorithm  can  be  solved  on  a  pro¬ 
grammable  scientific  calculator  and  would  provide  the  engineer  with  a  valuable 
predictive  tool  in  situations  where  he  or  she  has  no  access  to  a  computer  capa¬ 
ble  of  running  an  EOS-based  flash  routine. 
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Appendix  A 

MATHEMATICAL  PROOF  OF  THE  GENERALIZED  EQUATION 
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KNOWN: 


£  Z;  Ci  n  (1  +  a  Cj)  =  £  a<*-P>  £  z,-  C,  j  CN)  (A-l) 

*=i  p=i  «=i  l  j 

POSTULATE: 

£  zf  C,  fl  (1  +  a  Cj)  =  £  o<"-f>  {  n  ck  £  f(-ir'  £  777]  }a-2) 

i=l  y*i  p=  1  l*=l  >=1  L  ,=1  01  J  J 

=  £  a £  z.C,  j  ft  Ct  t  <-iy*‘  7^-  }  (A-3) 


k= 1  y=l 


Prove  that  the  coefficients  of  a  in  Equation  (A-l)  and  Equation  (A-3)  are 


equivalent: 


aN^cx,...,  c„..„  cA  =  ifi  Ck  £  (-iy+1 


PROOF: 


(A-4) 


We  can  express  the  a-function  as: 


N 

n  c* 

*=1 


&N-p(C Ci,...,  C^)  —  £  ap-i(C\ 1 Cj  »•••»  Ch  )  (A-5) 

N 

n  ch 

where  —p; —  represents  (N- 1)  terms:  N-p  =  (N- 1)  -  (p-1) 

W 

Eliminate  the  Cf1  term  in  the  RHS  of  Equation  (A-5)  by  rewriting  the  a-function 


ap-i(C f1 .....  Cf1 . C/71)  =  fy-ifCf1 C/71)  -  Cf1  ap_2(C f1 . C/7lj)(A-6) 
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ap-  2(^1 1 
ap-i(C  f 1 


C,-  C^1)  -  flp_2(C  j 1 -  C,-  1  Qp-?,(C  j 1 C^/ii  )(A-7) 

C,-  C/71)  =  flp-3(^i_1 »— »  C/71)  -  C,_1  flp_4(Cf] C^Ii  )(A-8) 


flp-(p-i)(C f1  *  Q  1»—»  Qv1)  -  1 1  Cn  !)  ~ 

cr^-"  a0(Cxl . C^ij)  (A-9) 

Recall  that  a0  =  1  and  then  substitute  Equations  (A-7),  (A-8),...t  (A-9)  into  Equation 
(A-6): 

Bp-iCCi1 Cj  Ctf1)  =  cip-\(Cil C^-1)  -  Ci*1  ap_2(Ci_1 Cyv1) 

+  C,  2  flp_3(C] 1 Cf 71)  -  C,  3  ap_4(C  1 1 C^-1)  + 

•••  ±  Ci-^-2)^^_l)(Cf,,...,Ci71)  ±  (A- 10) 

.  ..ter  writing  the  recursive  form  for  the  RHS  of  Equation  (A- 10),  the  equation 
becomes: 


ap.x(Cil,...,C-\...tCfil)  =  f  (-1  y*1  Ci-V-"  ap_j 

;=i 

Substitute  Equation  (A-ll)  into  Equation  (A-5): 


&N-p  1 » — »  Q  Ctf) 
Combine  C,  terms: 


N 

net 


k= 1 


C, 


t  (-iy+1 

>=1 


Vl 

r-l 


Q' 


aN_p(Ci„..,  CN)  =  n 

t=i 


t  (-ly*1 

j~\  Ci' 


(A-ll) 


(A- 12) 


(A- 13) 


Equation  (A— 13)  =  Equation  (A-4)  Q.E.D. 


Appendix  B 

ALGORITHM  FLOWCHART 
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Appendix  C 
COMPUTER  CODE 
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******************** 


7  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 


Advisor. 


M.S.  thesis 

Dr.  Michael  A.  Adewumi 


****++*+**++++*****+***++******+*****++*++*++++*+*+++++++++*++++++* 

*  Program  ALPHATEST  (FORTRAN  77) 
******************************************************************* 

* 

*  This  program  calculates  values  of  the  vapor  fraction, 

given  equilibrium  ratios,  Ki,  and  feed  mole  fractions, 
zi.  It  can  be  used  to  reproduce  experimental  results  of 
equilibrium  flashes. 

ALPHATEST  calls  ALPHACOEFF,  BUD  AN,  ALPHAPLOT,  and 

ALPHAROOT 

ALPHACOEFF  calls  subroutine  SYMFUNCTION 
SYMFUNCTION  calls  subroutine  DETERM  and  function  FACTOR 

VARIABLES:  alpha  =  calculated  system  vapor  fraction 

beta  =  experimental  system  liquid  fraction 
coefficient  =  coefficient  of  alpha  polynomial 
Ki  =  equilibrium  ratio  for  component  i 
molefrac  =  feed  mole  fraction  of  component  i 
Ncomp  =  number  of  components  in  feed 
Npress  =  number  of  data  sets  to  be  evaluated 
Pi  =  system  pressure,  psia 
Ti  =  system  temperature,  F 
xalpha  =  experimental  system  vapor  fraction 

It  is  formatted  to  input  zi,  temperature,  pressure,  liquid 
mole  fraction,  and  Ki 


IMPLICIT  REAL*8(a-h.o-z) 

REAL*8  Ki(500, 1 00) jaiolefrac(0: 1 00) 

PARAMETER  (Npress=  16  ,Ncomp=  1 0) 

DIMENSION  alpha(500),  beta(500),  coefficient(0:100), 
@  Pi(500),  Ti(500),  tarray(2),  xalpha(500) 


Data  Input 


The  number  of  components  (Ncomp)  and  the  number  of  data  sets 
to  be  run  (Npress)  are  specified  as  PARAMETERS' 


* 

* 


Open  and  Rewind  Input  and  Output  Files 


OPEN (uni  t=  1  ,file= 'indata ',status=  'old  0 
OPEN(unit=7,  file=  'table ',status=  'unknown ') 

OPEN(unit=8,file= 'plot ', statu  s= 'unknown  0 

RE  WIND(unit=  1 ) 

REWIND(unit=7) 

REWIND(unit=8) 

read(l,*)  (molefrac(i),  i  =  1,  Ncomp) 
do  1000  j  =  1,  Npress 

read(l,*)  Pi(j),  Ti(j),  beta(j) 
read(l,*)  (Ki(j,i),  i  =  1,  Ncomp) 
xalpha(j)  =  l.dO  -  beta(j) 

1000  continue 

* 

*  Choose  between  single  or  multiple  runs 

* 

write(6,*)  'Evaluate  one  data  set?  enter  1' 
write(6,*)  'Evaluate  all  data  sets?  enter  2' 
read(5,*)  numsets 

if(numsets  .EQ.  1)  then 

write(6,*)  'Enter  number  of  data  set  for  this  run' 
read(5,*)  j 
go  to  2100 
end  if 

do  2000  j  =  1,  Npress 

2100  write(7,*)  '  ' 
write(7  *)  '  ' 

write(7  *)  '  RUN  'j 

write(6  *)  'J  =  j 
write(7,2500)  Pi(j),Ti(j),beu(j) 

2500  form^( 'Pressure  =  ',f6.1,'  psia  Temperature  =  ',f6.1,'  F 
@Liquid  Mole  Fractiwi  =  ',f6.4) 


Call  subroutines 

Calculate  coefficients  of  polynomial 

call  ALPHACOEFF(NcompJ4pressj,molefrac,Ki,coefficient) 


♦ 

* 

♦ 


Predict  the  number  of  roots  on  [0,1]  by  Fourier-Budan  theorem 
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call  BUD AN(j>Jcomp,coefficientjium root) 


Solve  for  the  roots  by  Newton-Raphson  method 

call  ALPHAROOT(j,Ncomp,coefficient,xalphajiumroot, alpha) 

Generate  various  plots  (EDIT  the  file  to  remove  comments  for  specific 
options) 


call  ALPHAPLOT(NcompJjnolefrac,aJpha,coefficient,K/) 

■  ALPHAROOT  has  internal  output  section  to  compile  a  table 

'  listing  statistics  on  the  determination  of  alpha 

2000  continue 

******************************************************************** 

Produce  this  format  to  plot  data  points  as  dots: 

(PLOTFAT=20) 

2 

x(l)  y(l) 
x(l)  y(D 
2 

X(2)  y(2) 
x(2)  y(2) 
etc. 

******************************************************************** 
do  3000  j  =  1,  Npress 

write(8,3500)  alpha(j),xalpha(j),alpha(j),xalpha(j) 

3500  format( '2  V.el6.9,10x,el6.9y,e  1 6.9,10x.e  1 6.9) 

3000  continue 

CLOSE(unit=l) 

CLOSE(unit=7) 

CLOSE(unit=8) 

stop 

end 

******************************************************************** 


4  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
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*  University  Park,  Pennsylvania 

* 

*  M.S.  thesis 

* 

*  SUBROUTINE  ALPHACOEFF 

* 

*  This  subroutine  calculates  the  coefficient  for  each  tern  in 

*  the  general  polynomial  for  the  vapor  fraction,  alpha: 

* 

*  P(alpha)  =  cO  +  cl*alpha  +  c2*alpha**2  +  ...  + 

* 

*  c(Ncomp-l)*alpha*,*,(Ncomp-l) 

*  Equation  4.29  in  the  thesis. 

* 

********************* ******* ***************************************** 

SUBROUTINE  ALPHACOEFF(Ncomp(Npressjj>molefrac,Ki,coefficient) 

IMPLICIT  REAL*8(a-h,o-z) 

REAL* 8  Ki(500,100),  molefirac(0:100) 

INTEGER  p 

DIMENSION  coefficient^:  100),  c(100) 

OPEN(unit=  1 4,file=  'coeff' ,status=  'unknown  0 
OPEN(unit=  1 5 ,  file=  'coeff.plot  ',status=  'unknown  *) 

if(Ncomp  ,LT.  2)  then 

write(6,*)'You  cannot  flash  this  system' 
stop 
end  if 


Calculate  Ci  =  Ki  -  1 

do  0500  k  =  1,  Ncomp 
c(k)  =  Ki(jjjc)  -  l.dOO 
0500  continue 

p-loop  increments  the  power  of  alpha 


C  write(15,*)Ncomp 

do  1000  p  =  1,  Ncomp 
temporary  =  O.dOO 

do  2000  j  =  1,  p 

*  Zero-order  elementary  symmetric  function,  a0[l/Ci],  defined  as  1 

♦ 

if(p-j  .EQ.  0)  then 
apj  =  l.dOO 
go  to  2500 


u  u 


end  if 
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Call  subroutine  to  calculate  the  elementary  symmetric 
function,  apj 

call  SYMFUNCTION(Ncompj,p,c,apj) 

ratio  =  O.dOO 

do  3000  i  =  1,  Ncomp 

ratio  =  ratio  +  molefrac(i)/(c(i)**(j-l)) 
continue 

temporary  =  temporary  +  ((-l.dO)**(j+1))*apj*ratio 
continue 

coefficientCNcomp-p)  =  temporary 

write(14,*)'Coefficient(',Ncomp-p,')  =  ',coefficient(Ncomp-p) 
write(  1 5,*)Ncomp-p,coefficient(Ncomp--p) 

1000  continue 

return 

end 

******************************************************************* 

4  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 

M.S.  thesis 

SUBROUTINE  SYMFUNCTION 
This  subroutine  calculates  the  elementary  symmetric  function 

a(p— j){  1/Ci) 

******************************************************************* 


* 

* 

2500 

3000 

2000 


SUBROUTINE  SYMFUNCTION(Ncompj,p,c,apj) 

IMPLICIT  REAL*8(a-h,o-z) 

REAL*8  mmatrix(100,100) 

INTEGER  factor,p 


DIMENSION  c(100),  s(100) 
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*  Compute  the  power-sum  series:  s  =  sigma[  (1/Ci)**iambda  ] 


n  =  p  -  j 

do  1000  lambda  =  1,  n 
stun  =  O.dOO 

do  2000  i  =  1,  Ncomp 

sum  =  sum  +  (l.d0/c(i))**lambda 
2000  continue 

s(lambda)  =  sum 
1000  continue 

Build  the  matrix  MMATRIX 

do  3000  k  =  1,  n 

do  4000  1  =  1.  n 

if(l  .LE.  k)  mmatrix(kj)  =  s(k-l+l) 
if(l  .EQ.  k+1)  mmatrix(k,l)  =  DFLOAT(k) 
if(l  .GT.  k+1)  mmatrix(k,l)  =  O.dOO 
4000  continue 

3000  continue 

Since  al  { 1/Ci )  forms  a  [lxl]  matrix,  its  determinant  is  the 
element  itself 

if(p-j  .EQ.  1)  then 
det  =  mmatrix(l.l) 
go  to  5000 
end  if 

Compute  the  determinant  of  MMATRIX 


call  DETERM(mmatrix,n,det) 


Compute  the  elementary  symmetric  function 


5000  apj  =  det/factor(n) 

return 

end 

******************************************************************** 

* 

*  Function  to  compute  the  factorial 

* 

♦Jr****************************************************************** 
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FUNCTION  factor(n) 

INTEGER  factor.i.n 

factor  =  1 
if(n  .GT.  0)  then 
do  6000  i  =  241 
factor  =  factor*i 
6000  continue 
end  if 
end 

******************************************************************** 

* 

*  4  Dec  91 

* 

*  BRETT  D.  WEIGLE 

*  Petroleum  and  Natural  Gas  Engineering  Section 

*  Mineral  Engineering  Department 

*  College  of  Earth  and  Mineral  Sciences 

*  The  Pennsylvania  State  University 

*  University  Park,  Pennsylvania 

* 

*  M.S.  thesis 

* 

*  SUBROUTINE  DETERM 

* 

*  This  program  calculates  the  determinant  of  an  NxN  matrix. 

*  First,  partial  pivoting  is  performed  on  a  nonsingular  matrix  by 

*  Gaussian  elimination.  This  produces  a  triangular  matrix  whose 

*  determinant  can  be  calculated  by  computing  the  product  of  all 

*  the  diagonal  entries. 

*  The  augmented  matrix  does  not  contain  the  normal  last  column 

*  which  represents  the  right-hand  side  of  a  system  of  linear 

*  equations;  AUG  is  the  same  as  the  original  matrix. 

* 

*  VARIABLES: 

*  N  =  dimension  of  matrix 

*  AUG  =  augmented  matrix 

*  IJ,K  =  indices 

*  MULT  =  multiplier  used  to  eliminate  an  unknown 

*  PIVOT  =  used  to  find  nonzero  diagonal  entry 
******************************************************************** 

SUBROUTINE  DE'lERM ( augji.det) 

IMPLICIT  REAL*8(a-h,o-z) 

REAL*8  mult 
INTEGER  pivot 
DIMENSION  aug(100,100) 


* 

* 


Gaussian  elimination 


do  7000  i  =  1,  n 
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Locate  nonzero  entry 


if(aug(i,i)  .EQ.  0)  then 
pivot  =  0 
j  =  i  +  1 

3000  if((pivot  .EQ.  0)  .AND.  (j  .LE.  n))  then 

if(aug(j,i)  .NE.  0)  pivot  =  j 
j  =  j  +  1 
go  to  3000 
end  if 

if(pivot  .EQ.  0)  then 

print  'Matrix  is  singular' 
stop 
else 


Interchange  rows  I  and  PIVOT 


do  4000  j  =  i,  n 
temp  =  aug(ij) 
aug(i.j)  =  aug(pivotj) 
aug(pivotj)  =  temp 
4000  continue 

end  if 

end  if 


Eliminate  I-th  unknown  from  equations  1+1 .  N 


do  6000  j  =  i+1,  n 

mult  =  -aug(j,i)  /  aug(i.i) 

do  5000  k  =  i,  n 

aug(jjc)  =  aug(j,k)  +  mult  *  aug(ijc) 

5000  continue 

6000  continue 

7000  continue 

Calculate  the  determinant  of  matrix  AUG  by  computing  the 
product  of  the  diagonal  elements 

prod  =  l.dO 
do  8000  i  =  1,  n 
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do  9000  j  =  1,  n 

if(i  .EQ.  j)  prod  =  prod  *  aug(i,j) 
9000  continue 


8000  continue 


det  =  prod 


return 

end 

********************************************************************* 

* 

*  4  Dec  91 


BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 


* 

* 

* 

* 


M.S.  thesis 

SUBROUTINE  ALPHAROOT 


Subroutine  uses  an  interval-halving  technique  to  find 
the  best  root  value  to  initialize  the  Newton-Raphson  (N-R) 
iterative  calculations  which  determine  the  real  root  of 
the  alpha  polynomial  on  the  interval  [0.1]. 


* 

* 


PARAMETERS:  delta  =  alpha  increment 

epsilon  =  alpha  convergence  criterion 
VARIABLES:  alower  =  lower  bound  of  alpha  increment 
aupper  =  upper  bound  of  alpha  increment 
falpha  =  the  alpha  polynomial 
fprime  =  first  derivative  of  alpha  polynomial 
guess  =  iterative  variable  for  alpha 
guessO  =  initial  estimate  for  N-R 
intcount  =  #  of  intervals  until  sign  change 
iter  =  #  of  iterations  until  N-R  converged 
isign,isign2  =  flags  for  function  sign  change 
isign,isign2  =  flags  for  function  sign  change 
numroot  =  flag  for  #  of  zeros  (from  BUD  AN) 


********************************************************************* 


SUBROUTINE  ALPHAROOT(jNcomp,coefficient,xalpha, numroot, 
@  alpha) 

IMPLICIT  REAL*  8(a-h,o-z) 

INTEGER  p 

DIMENSION  alpha(500),  coefficient(0:100),  xalpha(500) 


PARAMETER(delta  =  O.OldO,  epsilon  =  l.d-06) 


Write  table  heading 


write(7,*)The  Fourier-Budan  Theorem  yields  \numroot,'  roots  on 
@this  interval' 
write(7,3500) 

3500  format('Intervals',4x,  Initial  Guess ',4x, 'Iterations ',4x, 'Calc. 

(3) Alpha', 4x, 'Exp.  Alpha  0 

Check  flag  NUMROOT  provided  by  subroutine  BUDAN  to  determine 
root-search  scheme 

if(numroot  .EQ.  0)  then 

write(6,*)  'No  root  on  the  interval  [0,1]  for  data  set  'J 
intcount  =  0 

write(7,3900)  intcount,xalpha(j) 

3900  format(i4,61x,£5.3) 

writeC?,*)^  root  on  the  interval  [0,1]' 
return 
end  if 

if(numroot  .EQ.  1)  then 
ilower  =  0 
iupper  =  0 
end  if 

if(numroot  .GE.  2)  then 
ilower  =  0 
iupper  =  1 
end  if 

Use  incremental  search  to  determine  initial  guess 
Interval  Endpoint  DO-Loop 

do  0400  jroot  =  ilower,  iupper 

intcount  =  0 

Test  the  polynomial  at  endpoint  for  initial  sign  value 


if(jrooL  EQ.  ilower)  then 
guess  =  DFLOAT(ilower) 
alower  =  guess 
aupper  =  alower  +  delta 
end  if 

if(jnx)L  EQ.  iupper)  then 
guess  =  DFLOAT(iupper) 


aupper  =  guess 
alower  =  aupper  -  delta 
end  if 
ichange  =  0 

falpha  =  O.dO 
do  1500  p  =  1,  Ncomp 

term  =  coefficient(Ncornp-p)*guess**(Ncornp-p) 
if(  (Ncomp-p)  .EQ.  0  )  term  =  coefficient^) 
falpha  =  falpha  +  term 

continue 

Initialize  ISIGN2  on  first  pass  with  endpoint 


if(ichange  .EQ.  0)  then 
if(falpha  .GE.  0.)  then 
isign2  =  1 
else 

isign2  =  0 
end  if 
end  if 


Note  the  sign  of  the  function 


if(falpha  .GE.  0.)  then 
isign  =  1 
else 

isign  =  0 
end  if 

Test  function  for  sign  change  and  increment  or  decrement  the 
search  variable  as  appropriate 

if(isign2  .EQ.  isign)  then 
if(jrcot  .EQ.  ilower)  then 
alower  =  aupper 
aupper  =  aupper  +  delta 
guess  =  aupper 
else  if(jroot  .EQ.  iupper)  then 
aupper  =  alower 
alower  =  aupper  -  delta 
guess  =  alower 
end  if 
end  if 


Exit  subroutine  if  no  sign  change  is  detected  on  interval  [0,1] 
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if(  (guess  .GT.  1.)  .OR.  (guess  .LT.  0.)  )  then 
write(6,*)  'No  root  on  the  interval  [0,1]' 
write(7,3800)  intcount,xalpha(j) 
format(i4,61x,f5.3) 

write(7,*)'No  root  on  the  interval  [0,1]' 
return 
end  if 


If  NO  sign  change  but  still  within  interval,  repeat  the  sequence 

if(isign  .EQ.  isign2)  then 
isign2  =  isign 
intcount  =  intcount  +  1 
ichange  =  1 
go  to  0600 

else 

If  there  IS  a  sign  change: 

Halve  the  interval  where  the  function  crosses  the  x  axis 

guessO  =  (alower  +  aupper)  /  2.d0 
end  if 


Provide  this  guess  to  Newton-Raphson  to  begin  calculations 


guess  =  guessO 


N-R  is  limited  to  1000  iterations  for  convergence 


iter  =  0 

do  1000  iteriimit  =  1,  1000 

iter  =  iter  +  1 
falpha  =  O.dOO 
fprime  =  O.dOO 

do  2000  p  =  1,  Ncomp 

falpha  =  falpha  +  coefficient(Ncomp-p) 

♦guess*  *(Ncomp-p) 

fprime  =  fprime  +  (Ncomp~p)*coefficient(Ncomp-p) 
*guess**(Ncomp-p-l ) 

continue 

calc  =  guess  -  falpha/fprime 
error  =  DABS((calc  -  guess)/calc) 
guess  =  calc 

if(error  .LE.  epsilon)  go  to  3000 


1000 


continue 


* 

* 

* 


print  ♦.'N-R  method  failed  to  converge  after  1000  iterations' 
Output  results  to  file  ’TABLE" 


3000  write(7,3600)  intcount.guessO, iter, guess, xalpha(j) 

3600  format(i4,13x,f5.3,10x,i4,13x,f9.6,7x,f5.3) 

alpha(j)  =  guess 

Begin  search  for  root  from  opposite  end  of  interval 


0400  continue 

return 

end 

*********************************************************************** 

6  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 

M.S.  thesis 

SUBROUTINE  ALPHAPLOT 

This  subroutine  is  used  for  several  purposes: 

1.  Plotting  F(alpha)  vs  alpha  [Rachford-Rice  obj  function] 

2.  Plotting  F(alpha)  vs  alpha  [polynomial] 

3.  Plotting  Fprime  vs  alpha  [polynomial] 

*********************************************************************** 


SUBROUTINE  ALPHAPLOTCNcompJ.molefrac, alpha, coefficient.Ki) 

IMPLICIT  REAL*  8  (a-h,o-z) 

REAL*  8  Ki(500, 1 00)  jnolefrac(  1 00) 

DIMENSION  alpha(500),coefficient(0: 1 00) 

INTEGER  p 

PARAMETER(start  =  O.OdO,  end  =  2.0d0,  stepsize  =  0.0005d0) 
OPEN(unit=l  1  ,ftfe='fa.  plot', status= 'unknown”) 


OPEN(unit=  1 2,file=  'fprime  .plot  ',status=  'unknown  0 


100 


Number  of  data  points  for  plotting 


number  =  IDINT((end  -  start  +  stepsize)/stepsize) 


************************************************************************ 

*  F(alpha)  vs  alpha  [polynomial] 

*  F'(alpha)  vs  alpha  [polynomial] 

*  Adjust  Ncomp,Npress  in  PARAMETER 

write(ll,*)  number 
do  1000  phase  =  start, end, stepsize 
falpha  =  O.dOO 
fprime  =  O.dOO 
do  2000  p  =  l.Ncomp 

falpha  =  falpha  +  coefficient(Ncomp-p)* 

@  phase*  *(Ncomp-p) 

C  fprime  =  fprime  +  ( Ncomp-p)*  coefficient comp-p )* 

C  @  phase**  (Ncomp-p-1) 

2000  continue 

write( 11,3600)  phase  .falpha 
C  write(llJ600)  phasefprime 

3600  format(f7.3.2x,f25. 12) 

1000  continue 

************************************************************************ 
C*  Rachford-Rice  objective  function 

C 

C  do  4500  k  -  Impress 

C  k  =  6 

C  write(ll,*)  number 

C  do  3000  phase  =  start.endstepsize 

C  falpha  =  O.dOO 

C  do  4000  i  =  Iff  comp 

C  falpha  =  falpha  +  (molefrac(i)*(Ki(k,i)  -  l.dO))  / 

C  @  (l.dOO  +  phase*(Ki(k,i)  -  l.dO)) 

C*  End  of  i  loop 

C  4000  continue 

C 

C  write(llJSOO)  phase  falpha 

C  3500  format(f732xf25.12) 

C*  End  of  phase  loop 

C  3000  continue 
C*  End  of  k  loop 


C  writef  11 3500)  phase  falpha 

C  3500  format(f73  2xf25 .12) 

C*  End  of  phase  loop 

C  3000  continue 
C*  End  of  k  loop 

C  4500  continue 

***************************************************** ****************** 


CLOSE(unit=ll) 

CLOSE(unit=12) 


return 

end 


********************************************************************* 

* 


*  6  Dec  91 

* 

* 

* 

* 

* 

* 

* 

* 

*  M.S.  thesis 

* 

*  SUBROUTINE  BUD  AN 

* 

*  Subroutine  uses  the  Fourier-Budan  Theorem  to  determine 

*  the  number  of  roots  that  the  alpha  polynomial  has  on  the 

*  interval  [u,v]. 

* 


BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 


*  PARAMETERS:  iu  =  lower  bound  of  alpha  interval 

*  iv  =  uppper  bound  of  alpha  interval 

*  VARIABLES:  coefficient  =  coefficient  of  alpha  polynomial 

*  dcoeff  =  coefficient  of  polynomial  derivatives 

*  deriv  =  derivatives  of  alpha  polynomial 

*  fvapor  =  the  alpha  polynomial 

*  ia,ib  *  #  of  sign  changes  for  derivative  series 

*  ivapor  -  alpha  =  vapor  fraction 

*  jsign.ksign  =  flags  for  derivative  sign  change 

*  numroot  =  number  of  zeros  on  the  interval 
********************************************************************* 


SUBROUTINE  BUDAN(J,Ncomp,coefficient,numroot) 

IMPLICIT  REAL*8(a-h,o-z) 

INTEGER  p 

DIMENSION  dcoeff(0: 100,0: 100),  coefficient(0:100),  deriv(0:100) 
PARAMETER(iu  =  0,  iv  =  1) 


C  DATA  ( coefficientfl ),  l  =  OF  comp-1)  /-l.fl.,-2.,3.,-4.J.  / 

OPEN(unit=2,file= 'test ',status=  'unknown  0 
REWIND(unit=2) 

ia  =  0 
ib  =  0 

do  0500  ivapor  =  iu,  iv,  1 

* 

*  Evaluate  the  polynomial  function  at  the  endpoints  iu  and  iv 


fvapor  =  O.dO 


do  0600  p  =  1,  Ncomp 

fvapor  =  fvapor  +  coefficient(Ncomp-p)*ivapor**(Ncomp-p) 
0600  continue 

write(2,*)  'fvapor  =  '.fvapor 
write(2,*)  ' 


Calculate  coefficients  of  first  derivative 


do  1000  n  =  Ncomp-1,  0,  -1 
dcoeff(0,n)  =  coefficient^) 
write  (2.*)  'dcoeff(0,',n,')  =  ',dcoeff(0ji) 

1000  continue 

write(2,*)  '  ' 

Calculate  coefficients  of  2nd-  and  higher-order  derivatives 
as  multiples  of  those  of  the  first  derivative 

do  1500  m  =  1,  Ncomp-1 

do  2000  n  =  Ncomp-m,  1,  -1 
dcoeff(m,n-l)  =  n*dcoeff(m-l,n) 
write  (2.*)  'dcoeff(\m,',',n-l,')  =  ', 

@dcoeff(mji-l) 

2000  continue 

write(2,*)  '  ' 

1500  continue 


Evaluate  the  derivative  series  at  the  endpoints  iu  and  iv 


do  3000  m  =  1,  Ncomp-1 
deriv(m)  =  O.dO 

do  4000  n  =  Ncomp-m,  1,  -1 

terra  =  dcoeff(m,n-l)*ivapor**(n-l) 
if(  (n-1)  .EQ.  0  )  term  =  dcoeff(m,n-l) 
deriv(m)  =  deriv(m)  +  tenn 
write(2,*)  'inter  derivf'.m,')  =  ',deriv(m) 
4000  continue 

write(2,*)  'total  deriv(',m,')  =  ',deriv(m) 
write(2,*)  ' 

3000  continue 


Count  the  sign  changes  between  the  terms  of  the  series 


ifffvapor.  LT.  0.)  then 
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ksign  =  0 
else 

ksign  =  1 
end  if 

write(2,*)  Tcsign  =  '.ksign/  for  fvapor' 

do  5000  i  =  1.  Ncomp-1 
if(deriv(i)  .LT.  0.)  then 
jsign  =  0 
else 

jsign  =  1 
end  if 

write(2,*)  'jsign  =  '.jsign,'  for  deriv(',i,')' 


* 

*  Increment  A  or  B.  depending  upon  the  endpoint  under  evaluation 

* 


if(ivapor  .EQ.  iu)  then 
if(ksign  .NE.  jsign)  then 
ia  =  ia  +  1 

write(2,*)  'ia  =  '.ia,'  for  deriv(',i,')' 
end  if 
end  if 

if(ivapor  .EQ.  iv)  then 
if(ksign  .NE.  jsign)  then 
ib  =  ib  +  1 

write(2,*)  'ib  =  '.ib,'  for  deriv(',i,')' 
end  if 
end  if 


ksign  =  jsign 

write(2,*)  Tcsign  =  '.ksign,'  after  deriv(',i,')' 
write(2,’*‘) 

5000  continue 
0500  continue 


* 

*  Pass  a  flag  to  calling  program  to  indicate  root  conditions 

* 


write(2,*)  'ia  =  '.ia,'  and  ib  =  '4b 

numroot  =  ia  -  ib 

write(2,*)  'numroot  =  '.numroot 

write(2,6000)  Ncomp-1,  numroot,  iu,  iv,  J 
6000  format(This  polynomial  of  order  '43/  has  '43,'  zeros  on  the  in 
@terval  ['42,742,']  for  J  =  '43) 

CLOSE(unit=2) 

return 

end 
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ABSTRACT 


The  Rachford-Rice  objective  function  for  flash  calculations  exhibits  a  nearly 
flat  slope  across  the  two-phase  region  and  sharp  discontinuities  near  the  dewpoint. 
These  features  make  iterative  solution  procedures  sensitive  to  the  initial  estimate  of 
the  root  and  prone  to  spurious  values  if  a  correction  step  throws  the  algorithm  out¬ 
side  the  two-phase  region  or  near  the  phase  boundary. 

This  work  centers  on  the  recasting  of  the  Rachford-Rice  objective  function  into 
a  polynomial  function  of  the  vapor  fraction,  a.  The  degree  of  this  polynomial  is 
one  less  than  the  number  of  components  in  the  system  and  its  coefficients  can  be 
calculated  from  the  feed  composition  and  .quilibrium  ratios.  A  recursive 
expression  is  developed  that  involves  symmetric  functions  and  can  be  easily  pro¬ 
grammed  on  a  computer  or  scientific  calculator. 

The  principal  advantage  of  this  new  form  of  the  objective  function  is  that  the 
theory  of  polynomials  is  well-developed.  The  location  of  their  zeros  can  be 
predicted  with  confidence  by  techniques  based  on  sound  mathematical  principles, 
such  as  the  Fourier-Budan  theorem.  The  a-polynomial  is  well-behaved  over  the 
two-phase  region  and  its  root  can  be  quickly  located  by  a  hybrid  method  of 
interval-halving  technique  and  Newton-Raphson  procedure.  The  validity  of  the  new 
objective  function  and  its  automatic  coefficient-generating  algorithm  are  tested  using 
several  multicomponent  systems  for  which  experimental  data  are  available. 

Results  of  these  tests  prove  conclusively  the  validity  of  the  generalized  polyno¬ 
mial  objective  function.  The  versatility  of  this  form  of  the  flash  objective  function, 
compared  with  the  original  Rachford-Rice  version,  is  demonstrated.  Another 


IV 


potential  advantage  of  the  polynomial  form  is  its  ability  to  handle  dilute  systems  in 
which  some  components  are  present  but  in  very  low  concentrations.  It  also  prom¬ 
ises  possible  usage  as  a  means  of  developing  appropriate  lumping  schemes. 
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NOMENCLATURE 


Roman 

A  =  pressure-dependent  constant  [Equation  (2.6)] 

A  =  parameter  for  the  PREOS  [Equation  (2.18)] 

A  =  number  of  sign  changes  in  derivative  series  (Section  5.2.4) 
a  =  intercept  of  log  Kp  vs  F  plot  [Equation  (2.3)] 
a  =  unitary  or  elementary  symmetric  function  [Equation  (3.41)] 
a(T)  =  attractive  constant  for  PREOS  [Equation  (2.11)] 
a{Tc)  =  attractive  constant  at  critical  point  for  PREOS  [Equation  (2.10)] 
B  -  parameter  for  the  PREOS  [Equation  (2.19)] 

B  -  number  of  sign  changes  in  derivative  series  (Section  5.2.4) 
b  =  translation  constant  for  Cox  chan  [Equation  (2.1)] 
b  =  molar  co-volume  for  PREOS  [Equation  (2.9)] 

C,  =  constant  used  in  the  objective  function  [Equation  (3.16)] 

C,  =  excluded  term  [Equation  (3.42)] 
c  =  slope  of  log  Kp  vs  F  plot  [Equation  (2.3)] 
c  =  volume  translation  parameter  [Equation  (2.21)] 

F  -  number  of  moles  in  feed  stream  (Section  3.1.1) 

F  =  component  characterization  factor  [Equation  (2.1)] 

Ki  -  equilibrium  ratio  (Section  1.2) 

L  =  number  of  moles  in  the  liquid  phase  (Section  3.1.1) 

M  =  interexpressibility  matrix  [Equation  (4.12)] 
m  =  element  of  matrix  M  [Equation  (4.31)] 

N  =  number  of  components  in  the  fluid  system  (Section  1.3) 


X 


n  =  number  of  terms  in  the  set  [Equation  (3.41)] 

P  =  a  partition  of  r  into  at  most  X  parts  [Equation  (4.8)] 

p  =  pressure 

Pi  =  pan  of  a  partition  [Equation  (4.1)] 

p(t)  =  polynomial  as  a  function  of  t  (Section  5.2.4) 

p(u)  =  polynomial  as  a  function  of  u  (Section  5.2.4) 

p  (v )  =  polynomial  as  a  function  of  v  (Section  5.2.4) 

R  =  real  gas  constant  [Equation  (2.8)] 

R  =  field  of  real  numbers  (Section  5.2.4) 
s  =  power  sum  symmetric  function  [Equation  (4.6)] 
sp  =  slope  of  plot  [Equation  (2.2)] 

T  =  temperature 

t  =  exponential  term  in  generating  function  [Equation  (3.41)] 
u  =  augmented  unitary  symmetric  function  [Equation  (4.7)] 
u  =  real  number  (Section  5.2.4) 

V  =  number  of  moles  in  the  vapor  phase  (Section  3.1.1) 

v  =  molar  volume 

v  =  pseudo  volume  [Equation  (2.21)] 

v  =  length  of  a  partition  [Equation  (4.2)] 
v  =  real  number  (Section  5.2.4) 
w  =  weight  of  a  partition  [Equation  (4.1)] 

X  =  fluid  "map"  coordinate  from  Varotsis  (Section  2.2) 

Xt  =  component  "map"  coordinate  from  Varotsis  (Section  2.2) 
x  =  mole  fraction  in  the  liquid  phase  (Section  3.1.1) 
x  =  argument  of  symmetric  function  [Equation  (3.41)] 

Y  =  fluid  "map"  coordinate  from  Varotsis  (Section  2.2) 
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y,  =  component  "map"  coordinate  from  Varotsis  (Section  2.2) 
y  =  mole  fraction  in  the  vapor  phase  (Section  3.1.1) 

Z  =  compressibility  factor  [Equations  (2.20),  (2.22)] 
z  =  mole  fraction  in  the  feed  stream  (Section  3.1.1) 


Greek 

a  =  vapor  fraction  [Equation  (1.1)] 

(3  =  isothermal  compressibility  (Section  2.3.2) 

5  =  binary  interaction  coefficient  [Equation  (2.16)] 

8  =  bulk  modulus  [Equation  (2.23)] 

t|  =  function  to  describe  a (T)  away  from  critical  point  [Equation  (2.1 1)] 

k  =  function  of  acentric  factor  in  PREOS  [Equation  (2.13)] 

Kg  =  function  of  acentric  factor  in  PRSV  EOS  [Equation  (2.26)] 

k,  =  parameter  in  PRSV  EOS  [Equation  (2.27)] 

p.  =  coefficient  of  polynomial  [Equation  (5.2)] 

7t,  =  multiplicity  of  a  part  in  a  partition  [Equation  (4.1)] 

V  =  EOS  variable  [Equations  (2.24),  (2.25)]] 

co,  =  Pitzer  acentric  factor  for  the  i-th  component  [Equation  (2.5)] 


Subscripts 

B  =  boiling  point  [Equation  (2.1)] 

c  =  critical  property 

i  J  JcJjn  =  individual  components  of  the  fluid  system 

ij  =  interaction  between  component  i  and  component  j  of  the  fluid  system 
k  =  convergence  pressure  (Section  2.2) 
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P  =  constant  pressure  [Equation  (4.29)] 

R  =  reduced  property 

RA  =  Rackett  compressibility  factor  [Equation  (2.22)] 

r  =  order  of  symmetric  function  [Equation  (3.41)] 

T  =  constant  temperature  [Section  2.3.2,  Equation  (4.29)] 

fluid  state  [Equations  (2.24),  (2.25)] 
component  of  fluid  system  [Equations  (2.24),  (2.25)] 
order  of  derivative 
first  derivative 
second  derivative 
third  derivative 

Abbreviations 

CPU  =  computer  central  processing  unit 

EOS  =  equation  of  state 

°F  =  degree  Fahrenheit 

LHS  =  left-hand  side  (of  an  equation) 

PREOS  =  Peng-Robinson  equation  of  state 
PRSV  EOS  =  Peng-Robinson-Stryjek-Vera  equation  of  state 
psia  =  pounds  per  square  inch,  absolute 
Q.E.D.  =  quod  erat  demonstrandum,  which  was  to  be  proved  (Appendix  A) 
RHS  =  right-hand  side  (of  an  equation) 

SRKEOS  =  Soave-Redlich-Kwong  equation  of  state 


SSM  =  successive  substitution  method 
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Chapter  1 


DISCUSSION  OF  THE  PROBLEM 


1.1  Introduction 


Determination  of  the  equilibrium  state  of  coexisting  liquid  and  vapor  phases, 
particularly  for  multicomponent  fluid  mixtures,  is  of  vital  interest  to  the  petroleum 
and  chemical  industries.  Many  processes  in  petroleum  production  and  refining 
involve  repetitive  flash  calculations  for  design  and  operational  purposes.  The  pri¬ 
mary  goal  of  performing  flash  calculations  is  to  determine  the  relative  amounts  and 
compositions  of  the  coexisting  phases  for  a  given  feed  composition  at  a  specified 
condition  of  temperature  and  pressure. 

This  work  is  confined  solely  to  two-phase  vapor-liquid  equilibrium  computa¬ 
tions,  although  its  results  will  no  doubt  find  application  in  multiphase  flash  prob¬ 
lems  in  the  future. 

1.2  The  Generic  Flash  Algorithm 

To  begin  the  calculation,  the  following  variables  must  be  specified:  the  system 
pressure  and  temperature,  the  molar  composition  of  the  feed  stream,  z, ,  and  an  ini- 

y. 

tial  estimate  of  the  equilibrium  ratios,  Af,  =  — .  The  process  is  assumed  to  occur 

**i 

under  isothermal  and  isobaric  conditions.  The  stages  of  the  calculation  are: 

1.  Compute  initial  estimates  of  the  equilibrium  ratios  by  one  of  the  established 
techniques  or  by  an  empirical  correlation. 
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2.  Calculate  the  phase  distribution  and  compositions  corresponding  to  the  given 
/sf-values.  This  involves  the  iterative  solution  of  the  following  objective  func¬ 
tion,  developed  by  Rachford  and  Rice  (1952): 


N 

L 


i=i 


Zj  (*,--!) 

1  +  a  (  K, ;  -  1  ) 


=  0 


where  a  is  the  vapor  molar  fraction. 


(1.1) 


3.  Use  an  equation  of  state  (EOS)  to  calculate  the  component  fugacities  in  each 
phase  and  check  for  equality. 

4.  If  equality  is  not  achieved  (i.e.,  the  phases  are  not  in  equilibrium),  correct  the 
AT-values  on  the  basis  of  the  fugacities  and  repeat  steps  2-4. 

This  correction  is  commonly  performed  using  a  successive  substitution-type 
method  or  a  second-order  Newton-type  scheme.  These  algorithms  are  well-known 
and  are  described  in  several  papers  [e.g.,  Risnes  et  al.  (1981);  Michelsen  (1982); 
Boston  and  Britt  (1978)]. 

Successful  implementation  of  the  generic  flash  algorithm  described  above 
requires  three  principal  elements.  These  are  (1)  a  general  estimate  of  the  set  of 
equilibrium  ratios  to  start  the  procedure;  (2)  a  good  equation  of  state  to  improve  Kt ; 
and  (3)  a  robust  objective  function  that  guarantees  convergence  to  a  single  value  of 
a.  A  poor  first  guess  of  Jf-values  may  produce  a  phase  split  that  is  physically 
impossible  under  the  prevailing  pressure  and  temperature.  Satisfactory  methods  are 
available  for  generating  these  values.  Furthermore,  existing  equations  of  state  do  a 
fairly  good  job  of  predicting  phase  properties,  and  other  efforts  continue  along  this 
line. 
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One  area  that  has  not  enjoyed  equal  amounts  of  attention  for  a  long  time  is  the 
form  of  the  objective  function.  Invariably,  the  Rachford-Rice  objective  function 
[Equation  (1.1)]  is  most  often  used.  Recent  investigations  (Warren,  1991)  have 
shown  that  this  objective  function  does  exhibit  some  strange  behavior  which  may 
affect  its  ability  to  generate  good  results  for  some  conditions. 

1.3  Objectives  of  the  Investigation 

This  study  is  aimed  at  evolving  a  generalization  of  the  new  polynomial  form 
of  the  Rachford-Rice  objective  function  developed  by  Warren  (1991).  The  resulting 
generalized  polynomial  function  of  the  vapor  fraction,  a,  should  be  capable  of  han¬ 
dling  an  ^-component  mixture.  The  coefficients  of  the  generalized  polynomial 
should  depend  on  only  two  variables,  the  molar  composition  of  the  feed  stream  and 
the  equilibrium  ratios,  and  should  be  easy  to  obtain,  either  analytically  or  numeri¬ 
cally.  Appropriate  algorithms  are  to  be  developed  for  this  purpose. 

The  principal  advantage  of  a  polynomial  form  of  the  flash-calculation  objective 
function  is  that  the  theory  of  polynomials  is  well-developed  and  semi-analytical 
solution  techniques  exist  for  equations  up  to  fifth-order  (Zaguskin,  1961).  For 
higher-order  polynomials,  the  Newton-Raphson  iterative  method  usually  provides  a 
fast  and  accurate  determination  of  the  roots. 

Determination  of  all  the  zeros  of  this  polynomial  is  unnecessary  since  the  phy¬ 
sics  of  the  problem  demands  that  only  the  zeros  on  the  bounded  interval  [0,1]  are  of 
practical  interest.  Furthermore,  the  physics  also  suggests  that  only  one  zero  (or 
value  of  a)  exists  on  this  interval,  which  represents  the  two-phase  vapor-liquid 
region.  It  can  be  shown  mathematically  that  this  is  indeed  the  case  for  well-defined 
systems,  as  will  be  demonstrated  in  5.2. 


Chapter  2 


LITERATURE  REVIEW 


A  survey  of  the  pertinent  literature  reveals  that  apparently  only  one  other 
worker,  Warren  (1991),  has  studied  the  particular  aspect  of  flash  calculations  tar¬ 
geted  in  this  research.  A  comprehensive  review  of  the  literature  pertaining  to  the 
use  of  cubic  equations  of  state  in  flash  calculations  was  conducted  in  order  to  pro¬ 
vide  a  reference  point  for  the  testing  of  the  polynomial  objective  function. 

This  review  is  sub-divided  into  three  sections:  flash  calculation  algorithms; 
equilibrium  ratios;  and  cubic  equations  of  state.  Particular  emphasis  is  laid  on  the 
Peng-Robinson  equation  of  state. 

2.1  Vapor-Liquid  Equilibrium  Flash  Calculations 

This  discussion  will  be  confined  to  two-phase  vapor-liquid  equilibria.  The 
work  to  date  concentrates  on  developing  robust  algorithms  with  rapid  convergence 
rates.  Robustness  implies  the  ability  to  continue  the  calculations  after  recovering 
from  a  spurious  value  of  the  vapor  fraction  computed  in  the  neighborhood  of  the 
critical  point  or  at  the  phase  boundaries.  Abhvani  and  Beaumont  (1987)  present  an 
excellent  review  of  EOS -based  flash  algorithms.  They  divide  the  papers  into  two 
categories  according  to  solution  method,  those  using  some  variant  of  the  successive 
substitution  method  (SSM)  or  those  employing  a  second-order  Newton-type  method. 

The  SSM  technique  is  the  traditional  solution  algorithm,  but  it  exhibits  a  poor 
rate  of  convergence  and  ’  stability  problems  close  to  saturation  points  and  in  the 
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critical  region.  Risnes  et  al.  (1981),  Michelsen  (1982),  and  Mehra  et  al.  (1983) 
made  attempts  at  acceleration  and  stabilization  of  this  method. 

Similarly,  many  workers  have  proposed  various  forms  of  second-order  Newton 
procedures  to  avoid  the  slow  rate  of  convergence  of  the  SSM,  such  as  Boston  and 
Brin  (1978),  Fussell  and  Yanosik  (1978),  Asselineau  et  al.  (1979),  Fussell  (1979), 
Baker  and  Luks  (1980),  and  Varotsis  et  al.  (1981).  Others  advocate  a  combination 
of  successive-substitution  and  Newton  methods;  the  former  is  used  to  provide  good 
initial  values  to  the  rapidly  converging  latter.  Informative  studies  include  Mott 
(1980,  1983),  Mehra  et  al.  (1982),  Michelsen  (1982),  Nghiem  et  al.  (1983),  and 
Abhvani  and  Beaumont  (1987). 

Benmekki  (1984)  developed  a  general  algorithm  for  flash  calculations  that  can 
utilize  any  cubic  equation  of  state  and  features  a  specified  calculational  path  for 
computing  the  phase  boundaries.  This  is  an  attempt  to  ensure  that  bubblepoint  and 
dewpoint  computations  originate  from  within  the  two-phase  region,  thus  guarantee¬ 
ing  meaningful  values  of  the  equilibrium  ratios. 

Warren  (1991)  made  a  radical  departure  from  previous  efforts  at  enhancing 
flash  calculation  algorithms  when  he  formulated  an  explicit  linear  equation  for  the 
vapor  fraction  of  a  binary  system.  He  successfully  extended  this  to  a  quadratic 
equation  for  a  ternary  system  and  a  cubic  equation  for  a  quaternary  mixture.  The 
success  achieved  by  Warren  and  the  possibility  of  the  existence  of  a  generalized 
polynomial  expression  for  the  vapor  fraction  in  a  two-phase,  N-component  fluid  sys¬ 
tem  motivated  the  current  work. 
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2.2  Vapor-Liquid  Equilibrium  Ratios 

The  use  of  initial  equilibrium  ratios  close  to  the  final  values  for  a  multicom¬ 
ponent  fluid  is  crucial  to  the  rapid  convergence  of  any  flash  calculation.  Experi¬ 
mental  values  are  preferred  because  the  prediction  of  Ki  for  a  particular  fluid  at 
various  combinations  of  temperature,  pressure  and  composition  requires  lengthy  cal¬ 
culations.  Therefore,  predictive  methods  for  A’-values  are  a  limiting  factor  in  the 
speed  and  robustness  of  any  flash  calculation  algorithm. 

The  expression  "equilibrium  constant"  was  coined  by  Souders  et  al.  (1932) 
and  was  defined  as  the  ratio  of  the  vapor  mole  fraction  to  that  of  the  liquid.  The 
basis  for  most  predictive  methods  had  its  genesis  when  Cox  (1923)  observed  that 
the  lines  on  a  semilogarithmic  plot  of  vapor  pressure  against  temperature  appeared 
to  converge  to  a  single  pressure.  Katz  and  Hachmuth  (1937)  demonstrated  an 
analogous  behavior  for  equilibrium  constants;  they  converged  to  unity  at  a  fluid 
mixture’s  critical  pressure. 

White  and  Brown  (1942)  attempted  to  develop  a  correlation  for  K-values  based 
on  this  "convergence"  pressure.  Hanson  and  Brown  (1945)  used  experimental  data 
to  correlate  the  convergence  pressure  (pk )  at  one  temperature  as  a  function  of  the 
molal  average  boiling  point  of  the  equilibrium  vapor  and  liquid.  They  also  showed 
that  the  convergence  pressure  concept  could  be  extended  from  binary  to  multicom¬ 
ponent  systems. 

Hadden  (1948,  1953)  produced  nomographs  for  equilibrium  constants  of  pure 
components  as  functions  of  temperature  and  pressure,  and  incorporated  convergence 
and  vapor  pressures  into  nomographs  for  mixtures.  He  demonstrated  that  mixture 
convergence  pressure  is  a  function  of  the  operating  temperature  and  of  the  liquid- 
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phase  composition  exclusive  of  the  lightest  component  concentration.  This  compo¬ 
sition  dependence  led  Muskat  (1949)  to  propose  the  use  of  the  term  "equilibrium 
ratio"  in  place  of  "equilibrium  constant."  Edmister  (1949)  presented  a  graph  involv¬ 
ing  the  ratio  of  differences  between  the  convergence  and  critical  pressures  and  the 
ratio  of  differences  between  the  system  and  critical  temperatures. 

Winn  (1952)  developed  nomographs  based  on  Hadden’s  (1948)  results  that 
allow  the  determination  of  /(-values  at  a  convergence  pressure  of  5000  psia.  For 
systems  with  pk  *  5000,  he  provides  a  translation  to  find  the  value  of  /(,  at  other 
"apparent"  convergence  pressures.  The  methods  proposed  by  these  three  authors 
require  chans  and  do  not  lend  themselves  to  computer  calculations. 

Hoffmann  et  al.  (1953)  attempted  to  extend  Cox’s  (1923)  vapor  pressure  graph 
for  the  purpose  of  determining  equilibrium  ratios  by  plotting  log  Kp  against  the 
component  characterization  factor  F,  where 


Kp  =  product  of  equilibrium  ratio  and  pressure 


b  =  constant  required  to  translate  the  vapor  pressure  curve 
for  a  hydrocarbon  onto  the  straight  line  of  the  Cox  chart 
TB  —  hydrocarbon  boiling  point 
T  =  system  temperature 


Brinkman  and  Sicking  (1960)  presented  an  iterative  method  for  finding  the  conver¬ 
gence  pressure  based  on  the  slope,  sP,  of  the  plot  mentioned  in  Hoffmann  et  al. 
(1953).  Then,  the  equilibrium  ratio  could  be  determined  as 


AT  =  —  c"F 
P 


(2.2) 
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Standing  (1979)  observed  that  the  composition  dependence  of  the  equilibrium 
ratio  is  negligible  at  pressures  below  1000  psia.  He  proceeded  to  combine  the  work 
of  Hoffmann  et  al.  (1953)  and  Brinkman  and  Sicking  (1960)  to  develop  a  correla¬ 
tion  for  AT-values  for  the  crude  oils  studied  by  Katz  and  Hachmuth  (1937): 


K  =  —  10(a  +  cF)  (2.3) 

P 

where  a  and  c  are  the  intercept  and  slope  (respectively)  of  log  Kp  vs.  F  plots  of  the 
abovementioned  oils.  Both  a  and  c  are  expressed  as  functions  of  pressure.  He  also 
presented  equations  for  the  heavy  fraction  and  the  common  reservoir  gases  N2,  C02 
and  H2S  (referred  to  as  permanent  gases). 

Wilson  (1969)  published  a  /f- value  equation  that  currently  enjoys  widespread 
use  in  flash  calculations: 

eB 

Kt  =  —  (2.4) 

PRi 

where 

B  =  5.37(1  +  (D,)(l  —  ~  )  (2.5) 

*Ri 

pRi  =  reduced  pressure  of  the  i  -th  component 
TRi  -  reduced  temperature  of  the  t  — th  component 
CD,  =  Pitzer  acentric  factor  of  the  i  -th  component 

Wilson’s  equation  fails  to  predict  accurate  equilibrium  ratios  for  most  fluids  above 
pressures  of  500  psia,  as  illustrated  by  Warren  (1991).  Whitson  and  Torp  (1981) 
attempted  to  correct  this  problem  by  re-introducing  the  system  convergence  pressure 
to  the  Wilson  equation: 


= 


Pci_ 

Pk 


A  -  1 


PRi 


(2.6) 
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where 


>4=1  — 


JL 

Pk 


14.7 


14.7 


0.6 


(2.7) 


pci  =  critical  pressure  of  the  i  -th  component 

Risnes  and  Dalen  (1984)  took  an  approach  based  on  the  equation  of  state  used 
in  the  flash  calculations.  Their  basic  idea  was  to  assume  the  mixture  or  feed  to  be 
liquid  and  then  evaporate  up  to  one-half  of  the  system  to  form  a  gas  phase  by  use 
of  the  fugacities.  The  initial  A'-values  then  could  be  calculated  from  the  resulting 
phases.  This  method  is  reported  to  perform  well  near  the  critical  point  and  along 
the  bubblepoint  line  but  often  fails  along  the  dewpoint  curve. 

Reportedly,  the  most  accurate  X-value  predictor  is  that  proposed  by  Varotsis 
(1989).  He  used  over  1000  experimental  equilibrium  ratios  to  construct  an  X-Y  plot 
such  that  each  reservoir  fluid’s  position  on  the  "map"  is  determined  by  its  coordi¬ 
nates  X  and  Y.  These  coordinates  are  described  by  a  polynomial  fitted  to  the 
apparent  pressure  mentioned  in  Winn  (1952).  He  proposes  an  equation  for  the  con¬ 
vergence  pressure  based  on  the  mole  fraction  of  methane  and  nitrogen  in  the  fluid. 

Each  pure  hydrocarbon  component  is  represented  on  the  map  by  its  own  set  of 
coordinates  (X, ,  K, ),  which  are  calculated  as  functions  of  the  component  acentric 
factor.  Specific  values  are  given  for  the  permanent  gases  and  correlations  based  on 
molecular  weight  are  specified  for  the  heptane-plus  fraction. 

Finally,  the  straight  line  that  joins  the  pressure  and  temperature  coordinates 
(X,  Y)  of  the  fluid  with  the  position  of  each  component  on  the  map  (X,  ,  X,  )  inter¬ 
sects  the  X-value  axis  at  a  point  that  corresponds  to  the  equilibrium  value  of  the 
selected  constituent.  Varotsis  (1989)  presents  tables  for  three  different  crude  oils 
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and  gas  condensates  at  varying  temperatures  and  pressures  that  display  AT-values 
remarkably  close  to  experimental  values.  They  are  an  order-of-magnitude  improve¬ 
ment  over  those  predicted  by  the  equations  of  Wilson  (1969)  or  Whitson  and  Torp 
(1981). 

The  method  of  Varotsis  (1989)  was  attempted  in  the  current  work.  His  K- 
value  predictor  was  formulated  using  data  for  crude  oils  and  gas  condensates  con¬ 
taining  the  C1  -  C6  alkane  series,  the  heptane-plus  pseudocomponent  and  the  per¬ 
manent  gases.  It  will  not  properly  describe  systems  (such  as  the  methane-ethane- 
propane  ternary)  containing  fewer  components  than  these  "typical"  reservoir  fluids. 
For  lack  of  a  suitable  replacement  expression  for  pk  ,  Wilson’s  equation  is  used  in 
the  current  work. 

2.3  Cubic  Equations  of  State  (EOS) 

The  equation  of  state  (EOS)  is  the  heart  of  a  modem  flash  calculation  algo¬ 
rithm.  Ideally,  it  should  be  able  to  accurately  represent  the  thermodynamic  proper¬ 
ties  of  the  fluid  of  interest  over  the  complete  range  of  operating  pressures  and  tem¬ 
peratures.  Since  engineering  applications  rarely  focus  on  an  isolated  chemical 
species,  the  EOS  should  incorporate  mixing  rules  that  allow  it  to  extend  its  predic¬ 
tive  capabilities  to  the  behavior  of  multicomponent  fluids.  Its  component-specific 
descriptive  parameters  should  be  readily  calculable  from  well-known  properties, 
such  as  critical  temperature  and  pressure,  molecular  weight  and  acentric  factor. 
Finally,  the  associated  computations  should  not  consume  excessive  computer  time, 
especially  if  the  equation  of  state  is  to  be  used  for  repetitive  calculations. 

The  engineer  is  faced  with  the  choice  of  using  a  complex  EOS  exhibiting  a 
high  degree  of  non-linearity  and  many  adjustable  parameters,  or  a  cubic  EOS  which 
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possesses  an  analytical  solution  and  requires  the  estimation  of  two  or  three  parame¬ 
ters.  Mathias  and  Benson  (1986)  presented  a  comparison  of  average  central¬ 
processing-unit  (CPU)  times  required  by  three  cubic  EOS  and  by  three  complicated 
EOS  to  compute  fugacity  coefficients  and  enthalpy  departures.  They  asserted  that 
the  time  required  for  any  of  the  candidate  equations  to  calculate  the  density  root  (or 
compressibility  factor  root)  is  negligible  compared  to  that  involved  in  executing  the 
various  mixing  rules.  In  fact,  for  systems  containing  more  than  about  six  com¬ 
ponents,  the  cubic  EOS  become  more  computationally  burdensome  than  the  compli¬ 
cated  EOS  simply  because  of  the  cross  terms  inherent  to  the  cubic  EOS  mixing 
rules. 

Engineers  frequently  use  cubic  EOS  because  they  work  well  over  the  range  of 
most  industrial  operating  conditions  and  are  easily  programmed  for  solution  on  a 
computer.  The  two  cubic  EOS  which  have  gained  the  widest  acceptance  are 
Soave’s  modifications  of  the  Redlich-Kwong  (1949)  equation  of  state  (SRKEOS) 
(Soave,  1972)  and  that  presented  by  Peng  and  Robinson  (1976b)  (PREOS).  The 
PREOS  and  suggested  improvements  are  examined  in  this  work  for  possible  use  in 
flash  calculations  because  of  the  author’s  familiarity  with  this  EOS. 

2.3.1  Development  of  the  Peng-Robinson  EOS 

Upon  the  success  of  the  SRKEOS,  Peng  et  al.  (1975)  undertook  a  further  study 
to  formulate  a  cubic  equation  of  state  with  an  improved  capability  to  predict  liquid 
densities  and  other  fluid  properties,  particularly  in  the  vicinity  of  the  critical  region. 

This  study  resulted  in  a  further  modification  of  the  attractive  pressure  term  of 
the  classical  equation  of  state  proposed  by  van  der  Waals  (1873).  The  result  was 
the  EOS  presented  by  Peng  and  Robinson  (1976b): 
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_  RT _ a(T) 

P  v-b  v(v+b)  +  biy-b) 

The  values  of  the  parameters  are  obtained  from 


(2.8) 


b  =  0.07780  — - 

(2.9) 

Pc 

R2Tc 2 

a(Tc)  =  0.45724  - 

(2.10) 

Pc 

a(T)  =  a (Tc)  •  r\(TR,  cu) 

(2.11) 

TV*  =  l+K(l-7>*) 

(2.12) 

k  =  0.3746  +  1.4850co  -  0.1644C02  +  0.01667a)3 

(2.13) 

Equation  (2.12)  has  the  same  form  as  that  used  by  Soave  (1972)  but  k  was  obtained 
by  fitting  a  larger  range  of  vapor  pressure  data  as  a  function  of  the  reduced  tem¬ 
perature  and  the  acentric  factor  (Pitzer  et  al.,  1955)  of  each  component. 

In  order  to  use  the  equation  for  systems  containing  more  than  one  component, 
the  following  mixing  rules  are  presented: 


a  =YLxixjaij  (2.14) 

*  j 

b-'Zxibi  (2.15) 

I 

where 

««;  =0  ~  (2.16) 

Equations  (2.14)  and  (2.15)  are  consequences  of  the  mixing  rule  proposed  by  Kay 
(1936),  while  Equation  (2.16)  was  developed  by  Zudkevitch  and  Joffe  (1970).  The 
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experimentally  determined  binary  interaction  coefficient,  81; ,  characterizes  the 
binary  formed  by  component  i  and  component  j.  The  importance  of  8i;  in  accu¬ 
rately  reproducing  P-V-T  data  was  discussed  by  Peng  and  Robinson  (1976b)  and  by 
Robinson  et  al.  (1985). 


The  PREOS  can  be  written  in  the  form  of  a  cubic  equation  in  the  compressibility 
factor 

Z3  -  (1  -  B  )Z2  +  (A  -  3B2  -  2B  )Z  -  {AB  -  B2  -  B3)  =  0  (2.17) 


where 


A  = 


r2T2 


B  —  bp* 

RT 


Z  =  ^ 
RT 


(2.18) 

(2.19) 

(2.20) 


2-3-2  Selection  of  the  Proper  Root  in  Cubic  EOS 

Equation  (2.17)  yields  one  or  three  roots  depending  upon  the  number  of  phases 
in  the  system.  The  authors  stated  that,  in  the  two-phase  region,  the  largest  root  is 
for  the  compressibility  factor  of  the  vapor  while  the  smallest  positive  root 
corresponds  to  that  of  the  liquid. 

Lawal  (1987),  however,  asserted  that  this  criterion  was  insufficient  to  select  the 
proper  root.  He  proved  that,  in  the  event  of  multiple  real  roots,  the  smallest  of  the 
positive  roots  larger  than  or  equal  to  B  must  be  chosen  for  the  compressibility  of 
the  liquid. 
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Asselineau  et  al.  (1979)  compared  the  calculated  volume  to  the  pseudo-critical 
volume  to  assign  the  root  to  the  proper  phase,  under  specific  conditions.  Poling  et 
al.  (1981)  examined  the  order  of  magnitude  of  the  isothermal  compressibility, 
(3  =  -(dv  /dp  )j  /v ,  to  ascertain  the  presence  of  the  liquid  or  vapor  phase.  Gosset  et 
al.  (1986)  offered  two  discriminants,  one  based  on  the  Cardan  criterion  for  the 
number  of  real  roots  for  a  cubic  equation  and  a  heuristic  approach  similar  to  that  of 
Asselineau  et  al.  (1979). 

2.3.3  Modifications  to  the  Peng-Robinson  EOS 

Numerous  attempts  have  been  made  to  correct  for  the  deficiencies  inherent  in  a 
cubic  equation  of  state  by  introducing  additional  parameters  into  the  PREOS. 
These  changes  improve  some  aspect  of  the  EOS’s  performance  (usually  liquid  den¬ 
sity  predictions)  but  at  the  cost  of  increased  complexity  and  the  requirement  for 
tables  or  correlations  to  determine  the  additional  parameter(s)  for  each  fluid  com¬ 
ponent.  This  review  will  touch  on  a  limited  number  of  these  studies. 

2J.3.1  Volume  Corrections 

The  modification  of  the  SRKEOS  proposed  by  PSneloux  et  al.  (1982)  also 
formed  the  basis  for  two  other  studies  concerned  with  the  PREOS.  These  authors 
suggested  that  the  use  of  a  "pseudo  volume"  defined  by 

v  =  v  +  £  cixi  (2-21) 

< 

could  be  used  to  effect  a  translation  along  the  volume  axis,  leaving  unchanged  the 
predicted  equilibrium  conditions.  They  chose  c  so  that  correct  saturated  liquid 
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densities  were  exactly  reproduced  at  the  reduced  temperature  TR  =0.7.  They 
rejected  the  acentric  factor  as  a  correlating  parameter  in  favor  of  the  Rackett 
compressibility  factor,  Zra  ,  developed  by  Spencer  and  Danner  (1972): 


c  =  0.40768 — —(0.29441  -  ZM )  (2.22) 

Pc 

Their  third  parameter  did  improve  predictions  of  saturated  liquid  densities. 

Almost  simultaneously,  Jhaveri  and  Youngren  (1988)  and  Mathias  et  al. 
(1989)  presented  three-parameter  modifications  of  the  PREOS  based  on  the  work  of 
Peneloux  et  al.  (1982).  The  first  authors  correlated  the  third  parameter,  c,  with 
molecular  weight.  The  second  study  retained  the  Peneloux-Rauzy-Freze  volume 
correction  scheme  but  added  a  further  term  involving  the  bulk  modulus  to  handle 
the  critical  region.  The  bulk  modulus  is  dimensionless  and  is  defined  as: 


5 


(2.23) 


From  an  examination  of  the  graphs  accompanying  both  publications,  the  work  of 
Mathias  et  al.  (1989)  seems  to  produce  results  closer  to  the  experimental  values  for 
saturated  volumes  and  densities. 


23.3.2  Temperature  Dependence 

Xu  and  Sandler  ( 1987a, b)  postulated  that  the  molar  co- volume  term,  b,  is  not 
independent  of  temperature  and  they  disputed  the  fitting  of  vapor  pressures  used  by 
Peng  and  Robinson  (1976b)  to  characterize  the  attractive  constant,  a.  They  corre¬ 
lated  vapor  pressure  and  volume  data  for  16  components  at  both  subcritical  and 
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supercritical  conditions  and  proposed  to  replace  the  numeric  coefficients  of  a  and  b 
found  in  Equation  (2.9)  and  Equation  (2.10)  with: 

(2.24) 

1=0 

and 

Vi1  =  (2.25) 

i=0 

where  i  refers  to  the  species  and  I  denotes  either  subcritical  or  supercritical  condi¬ 
tions. 

Wu  and  Sandler  (1989)  generalized  the  temperature-dependent  parameters  of 
Xu  and  Sandler  (1987a,b)  by  performing  least-squares  fits  of  and  as  func¬ 
tions  of  acentric  factor  and  reduced  temperature.  They  were  able  to  accomplish  this 
task  only  for  the  n-alkane  series  because  of  insufficient  data.  For  their  intended 
application  of  the  work  (petroleum  reservoir  simulation),  they  envisioned  the  use  of 
the  fluid-specific  parameters  for  the  permanent  gases,  water  and  light  ends  and  the 
generalized  parameters  for  the  heavy  pseudocomponents. 

Stryjek  and  Vera  (1986a, b,c,d)  re-worked  Equation  (2.13)  to  obtain  a  better 
reproduction  of  vapor  pressure  data  at  low  reduced  temperatures: 

Kq  =  0.378893  +  1.48971530)  -  0. 1 7 1 3 1 848co2  +  0.0196554©3  (2.26) 

and  modified  Equation  (2.12)  by  the  introduction  of  one  compound-characteristic 
adjustable  parameter,  Kj: 

k  =  Kq  +  Kj  (1  +Tr  *)(0.7  -  7> )  (2.27) 

Stryjek  and  Vera  (1986b)  and  Proust  and  Vera  (1989)  listed  values  of  Kj  for  over 
160  compounds  of  industrial  interest.  Stryjek  and  Vera  (1986d)  and  Wilczek-Vera 
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and  Vera  (1987)  examined  mixing  rules  of  varying  complexity  foi  u«r  v,ith  the 
Peng-Robinson-Stryjek-Vera  (PRSV)  EOS.  For  the  current  work,  the  PRSV  EOS 
with  the  original  PREOS  mixing  rules  (as  formulated  by  Stryjek  and  Vera,  1986b) 
is  used  and  produces  noticeably  better  results  than  did  the  PREOS. 

2.3.4  References  on  Cubic  EOS 

Abbott  (1979)  and  Martin  (1979)  presented  comprehensive  reviews  of  cubic 
equations  of  state  available  at  that  time,  and  Vidal  (1983)  and  Vera  et  al.  (1984) 
updated  the  topic.  Huron  and  Vidal  (1979)  proposed  composition-dependent  mixing 
rules  while  Mathias  and  Copeman  (1983)  discussed  mixing  rules  dependent  on 
volume.  Finally,  Peng  and  Robinson  (1976b),  Peng  and  Robinson  (1977),  Robinson 
and  Peng  (1978),  Robinson  (1979)  and  Peng  (1986)  developed  specific  applications 
of  their  EOS. 


Chapter  3 


DEVELOPMENT  OF  THE  POLYNOMIAL  FUNCTION 
FOR  SIMPLE  SYSTEMS 


This  chapter  discusses  the  work  published  by  Rachford  and  Rice  (1952)  and 
Warren  (1991)  on  performing  flash  calculations.  It  shows  the  development  of  the 
Rachford-Rice  objective  function  [Equation  (1.1)]  and  extends  Warren’s  work  as  a 
precursor  to  developing  a  generalized,  multicomponent  equation  for  the  vapor  frac¬ 
tion. 

3.1  The  Rachford-Rice  Flash  Objective  Function 

We  will  briefly  examine  the  derivation  of  the  Rachford-Rice  objective  function 
that  is  universally  used  today  in  flash  calculations.  After  plotting  its  behavior,  it 
will  become  plain  why  it  is  so  difficult  to  solve  by  iterative  techniques  such  as  the 
Newton-Raphson  method. 

3.1.1  The  Material  Balance  Development 

Flash  calculations  are  used  to  determine  the  compositions  and  quantities  of  the 
vapor  and  liquid  phases  at  equilibrium  which  result  when  an  N-component  fluid  of  a 
particular  composition  is  subjected  to  a  particular  pressure  and  temperature.  The 
composition  of  the  feed  stream,  F,  is  denoted  by  1  z,  and  it  flashes  into  L  moles  of 
liquid  with  composition  Z  x, ,  and  V  moles  of  vapor  with  composition  I  y(  .  The 
resulting  material  balance  equations  are: 
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F  =  L  +  V 

(3.1) 

Fzt  =  Lxt  +  Vyi 

(3.2) 

As  defined  in  Chapter  1,  the  equilibrium  ratio  is: 


Ki  = 


yj_ 

x> 


(3.3) 


and,  rearranging,  one  is  left  with  the  equation: 


y,  =xlKi 


(3.4) 


Substituting  Equation  (3.4)  into  Equation  (3.2)  yields: 


Fzi  =  Vxi  Ki  +  Lxi 


(3.5) 


Simplify  by  isolating  the  x,-  term  and  dividing  through  by  F: 


Zi  =  X; 


(3.6) 


Dividing  Equation  (3.1)  through  by  F  and  solving  for  —  yields: 


(3.7) 
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Substituting  Equation  (3.7)  in  Equation  (3.6)  and  simplifying  the  equation  results  in: 


Imposing  the  constraint  of  £  jc,  =  1  on  Equation  (3.8)  leaves: 


«=i  !  + 


Rearranging: 


-1=0 


1  +  j  (  Ki  -  1  ) 


Recalling  Equation  (3.4),  we  can  write: 


»  = 


1  +  j  (  AT,  -  1  ) 


Imposing  the  constraint  of  £  y,-  =  1  on  Equation  (3.11)  yields: 


(3.10) 


(3.11) 


-1=0 


1  +  J  (  Ki  ~  1  ) 


(3.12) 
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Combining  Equation  (3.10)  and  Equation  (3.12)  leaves: 


N 

2 


i=l 


Zi  (  Kj  ~  1  ) 

1  +  -£  ( AT,  -  1  ) 


=  0 


Defining  the  vapor  fraction,  a,  as: 


(3.13) 


(3.14) 


and  substituting  Equation  (3.14)  into  Equation  (3.13)  yield  the  Rachford-Rice  objec¬ 
tive  function: 


N 

2 


/=! 


2,-  (  Kj  -  1  ) 

1  +  a  (  Ki  -  1  ) 


(3.15) 


3.1.2  A  Graphic  Representation  of  the  Rachford-Rice  Objective  Function 

As  the  authors  noted,  their  formulation  of  the  objective  function  was  prone  to 
instability  near  the  values  of  a  that  represented  the  phase  boundaries,  namely,  0  and 
1.  They  showed  that  the  slope  of  the  function  near  these  points  may  be  quite  steep. 
It  is  this  feature  that  tends  to  throw  derivative- based  root-finding  techniques  out  of 
the  two-phase  region,  yielding  spurious  roots. 

Figure  3.1  depicts  the  behavior  of  the  objective  function  over  a  wide  range  of 
a  for  a  binary  system  of  70%  methane  and  30%  ethane  (Bloomer  et  al.,  1953). 
Figure  3.2  does  the  same  for  a  ternary  system  consisting  of  85%  methane,  10% 
ethane  and  5%  propane  (Parikh  et  al.,  1984).  Although  values  of  the  vapor  fraction 
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have  no  physical  meaning  outside  the  interval  [0,1],  these  graphs  serve  to  illustrate 
how  ill-behaved  the  objective  function  is.  Its  slope  is  almost  flat  as  it  traverses  the 
two-phase  region  and  it  is  plagued  by  spiky  singularities. 

This  work  will  attempt  to  develop  a  new  expression  for  a,  one  that  possesses 
reasonable  slope  over  the  desired  interval  and  has  no  discontinuities  near  the  phase 
boundaries. 

3.2  Warren’s  Explicit  Equations  for  the  Vapor  Fraction 

Warren  (1991)  expanded  the  Rachford-Rice  objective  function  into  a  polyno¬ 
mial  in  a  for  a  binary,  ternary  and  quaternary  fluid  system.  He  did  this  by  setting 
N  equal  to  2,  3  or  4,  respectively,  and  reducing  the  resulting  equations  to  their  sim¬ 
plest  polynomial  form  by  algebraic  manipulations.  To  demonstrate  the  validity  of 
his  work,  Warren  also  showed  that  the  higher-order  polynomials  would  reduce  to 
those  for  smaller  systems  when  the  appropriate  mole  fractions  and  equilibrium  con¬ 
stants  were  removed. 

We  will  assume  (as  did  Warren)  that,  under  isobaric  and  isothermal  conditions, 
the  equilibrium  constant  does  not  change  such  that  the  quantity  (K,  -  1 ),  which 
appears  in  the  objective  function,  can  be  represented  by  a  constant,  C; . 

We  will  reproduce  the  entire  process  here  for  a  binary  system  but  will  show 
only  the  final  result  for  a  ternary  and  quaternary  system,  since  the  algebra  can  be 
quite  tedious  and  repetitive. 
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3.2.1  Binary  System 


Starting  with  the  objective  function  as  defined  in  Equation  (3.15): 


£  . ,o 

itt  1  +  a  ( K, |  -  1  ) 


(3.15) 


and  defining  C,  =  AT,  -  1 ,  Equation  (3.15)  can  be  rewritten  as: 


N  Z ;  C, 


y  — - — - —  =  o 


(3.16) 


For  a  two-component  system,  setting  N  =  2  in  Equation  (3.16)  and  expanding 
term-wise  yields: 


- i — 1—  +  — t— t—  =  0 

1  +  aC,  1  +  a  C2 


(3.17) 


Moving  the  terms  with  the  subscript  ”2"  to  the  right-hand  side  of  the  equation: 


1  +  a  Cx 


z2  ^2 
1  +  a  C2 


(3.18) 


By  multiplying  each  side  by  (1  +a  Cj)  (1  -kx  ,  one  obtains: 


(rj  C,)  (1  +  a  C2)  =  -  (i 2  Cj)  (1  +  a  Cx) 


(3.19) 


Expanding  each  side  yields: 


Zj  Cj  +  CXZj  Cj  C  2  =  ~~  2  2  C  2  ~  Oi  Z2  C  j  C  2 


(3.20) 
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Move  a  terms  to  the  left-hand  side  of  the  equation  and  all  remaining  terms  to  the 
right-hand  side,  then  recall  that,  for  a  binary  system,  Zj  +  z2  =  1 : 


a(C,C2)--(x1C1  +  r2C2)  (3.21) 

Dividing  both  sides  through  by  Cj  C2  and  substituting  (A',-  -  1)  yields  the  explicit 
form  of  the  objective  function  for  a  binary  system: 


*i  zi 

K2  -  1  +  Kx  -  1 


(3.22) 


3.2.2  Ternary  System 


or  +  a 


3 

X  2« 

1=1 


r  _  *\ 


X<v 


iH 

j*> 


3 

+  X 

1=1 


=  0 


(3.23) 


3.2.3  Quaternary  System 


a3  +  a2 


4  (1  - 

h  ci 


+  a 


4 

X*i 

i=l 


X  cj 

L /*' 


nc; 


4  zi 

+  It1- 

i=1  n  cy 

j* 


=  0 


(3.24) 


3  J  Extension  of  Warren’s  Work  to  Larger  Systems 


Warren’s  method  can  be  used  to  develop  polynomial  expressions  for  systems 
having  five,  six  and  seven  components.  It  will  be  observed  that  the  terms  of  the 
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equation  expand  in  a  regular  fashion,  thereby  suggesting  the  possibility  of  develop¬ 
ing  a  recursive  relationship  dependent  only  on  N,  the  number  of  components.  Only 
the  algebra  for  the  five-component  system  will  be  presented,  as  those  for  six-  and 
seven-component  systems  follow  the  same  procedure. 


ZiCi  z2^2  z3  ^3  z4  ^4  z5  ^5 

1+aC1!  l+aC2  l  +  aC3  1  +  a  C4  1  +  a  C5 

5 

Multiplying  through  by  I"[  (1  +  a  C,  )  yields: 

i= l 


(3.16) 


(3.25) 


z i  Cj(  1  +  oc  C2)  (1  +  ot  C 3)  (1  +  cc  C 4)  (1  +  ot  C 5)  + 

z2  C2  (1  +  a  Cj)  (1  +  a  C3)  (1  +  a  C4)  (1  +  a  C5)  + 

z3  C3  (1  +  a  C])  (1  +  a  C2)  (1  +  a  C4)  (1  +  a  C5)  + 

z4  C4  (1  +  a  Cj)  (1  +  a  C2)  (1  +  a  C3)  (1  +  a  C5)  + 

z5  C5  (1  +  a  Cj)  (1  +  a  C2)  (1  +  a  C3)  (1  +  a  C4)  =  0  (3.26) 

Expanding  each  term: 

a4Zj  Ci  C2C3  C4C5  + 


o?ZiCi(C2C^C^  +  C2C2C^  +  C2C^C^  +  C2C^C^)  + 
cx2  Z\  Ci  (C2C3  +  C2C4+C2C5  +  C3C4+  C3C5  +  C4C5 )  + 
cc  z  1  Ci  (C2  +  C3  +  C4  +  C5  )  +  z  1  Ci  + 
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a4  Z4C1  C2C3C4C5  + 

a3  z  4  C  4  (C 1  C2C3  +  C1  C2C5  +  C2C3C5  +  C1  C  3  C  5  )  + 
a2  z4  C4  (Ci  C2  +  C1  C3  +  C1  C5  +  C2C3  +  C2C5  +  C3C5)  + 

OtZ4C4(Ci  +  C2  +  C3  +  C5)  +  Z4C4  + 


CX4  Z5C1C2C3C4C5  + 

CX3  Z  5  C  5  (C 1  C2C3  +  C1  C2C4  +  C1  C3  C4  +  C2C3C4)  + 

(X2  Z  5  C  5  (C 1C2  +  C1C3  +  C1C4  +  C2C3+  C2C4  +  C3C4)  + 
(X  Z  5  C  5  (C  i  +  C2"t"C3  +  C4)  +  Z5Cj  —  0 


(3.27) 


5 

Dividing  through  by  the  term  J~[  C,  and  adding  like  terms  yields: 

i=l 
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a2« 


C  ]  C2  C3  C4 


1  +^  +  ^4-  +  ^  +  _i^+  1 


C2  C3  C 2  C4  C2  C5  C3  C4  C3  C5  C4  C« 


1  +^  +  _V  +  _i_+  1 


C 1  C3  Cj  C4  C 1  C5  C3C4  C3C5  C4C5 


1  +-^+^+_I_+^.+  1 


CiC2  Cj  C4  Ci  c5  c2c4  c2  c5  c4c« 


1  +^L  +  ^  +  _l_  +  -L-+  1 


Ci  C2  C1C3  Ci  C5  C2C3  c2c5  C3  Cj 


1  +^Lr.  +  ^_  +  ^+  >  +  1 


C1C2  Ci  C3  CtC4  C2  C3  C2  C4  C3C4 


ai  z, 


■  C2  +  C3  +  C4  +  C5  ■ 

x  7 

'  Cj  +  C3  +  C4  +  C5  ' 

C 2  C3  C4  C 5 

t  z  2 

Ct  c3  c4  c5 

C 1  +  C2  +  C4  +  C5 
Ci  C 2  C4  C5 

C1  +  C2  +  C3  +  C4 

Ci  C2  C3  C4 


+  z4 


C 1  +  C2+C3  +  C5 

Ci  c2  c3  c5 


t  *1  (  *2 
C2C3C4C5  c  1  C3  C4  C5 


23  z4  z5 

c,c2c4  c5  +  c,c2  c3  c,  +  c,c2c3c4  =  0  (3-28) 

To  maintain  similarity  with  the  forms  of  the  quaternary  and  ternary  equations,  we 
can  separate  the  general  term  in  the  coefficient  for  a  in  Equation  (3.28)  into  four 
fractions: 


C;  +  C*  +  C/  +  cm 

1 

1 

1 

1 

Cy  C*  C/  Cm 

-  Z; 

C*  C/  Cm 

L 

Cj  Cl  Cm 

Cy  C*  Cm 

| 

Cy  C*  C/ 

(3.29) 
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After  multiplying  each  of  the  five  fractions  by  z,  and  collecting  terms  with  common 
denominators,  the  following  form  appears: 


Z; 


*■ +  zi 


Cl  Cl  C„ 


(3.30) 


5 

We  can  invoke  z,-  =1  to  construct  fractions  with  similar  subscripted  terms  in 

i= l 


both  numerator  and  denominator: 


1-Zk-Zi  -zm 
Cfc  C/  Cm 


(3.31) 


This  yields  a  final  polynomial  expression  of  the  same  general  form  as  those  of  War 
ren  (1991): 


a4  +  a3  i 


1  ~zi 

Ci 


+ 


1  ~  z2 
C2 


+ 


1  -  z3 

C3 


Izii  +  iziil 

c4  c5 


1  —  z  1  —  z2  1-zl_z3  1~Z1-Z4  1  ~  z  1  —  z5  1  —  z2  —  z3 

C\C2  C\  C3  Cj  C4  ^iC5  C2C3 


1  —  z2  —  z4  1  —  z2  —  z5  1  z  3  z  4 

C 2  C 4  C 2  C 5  C3C4 


*  ~  z3  ~  z5 

c3c5  + 


1  -  z4  ~z5 

c4  c5 


a 


1  ~  Z  1  ~  Z2  -  z3  1  -  Zt  -  Z2  -  Z4  1  -Z;  -Z2-Z5  1  ~Zt  ~Z3~Z4 

CjC2C3  CjC2C4  CiC2C5  C,C3C4 


1  -Zi  -Z3-Z5  1  -  Zj  -  z4  -  z5  l-z2-z3-z4  1  -  z2  -  z3  -  z5 

+ - — - + - — — — - + - -  —  - +  • 


C \C  2C  $  CiC,C5 

l-z2-z4-z5  l-z3-z4-z5 


C  2C  4C  5 


C3C4C5 


c  2C  3C  4 
Z1 

C2C3C4C5 


C2C3C5 
z2 

Cj  c3  c4  c5 


C1C2C4C5  C>  C2C3  C5  C!C2C3  c4 


=  0 


(3.32) 
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This  yields  the  polynomial  expression  for  a  quinary  system: 


a4  +  a3 


N* 

l 

i 

+  a2 

4  5  (1  -Zi-Zj) 

1/5  c,  J 

T  U 

[m  ;=7+ 1  C<  Cy  J 

fit  ^cVJk) 

L«=l  ;'=»+!  k=i+2  W  *■'* 


5 

+  5  5 

,=1  n  C) 

j** 


=  0 


(3.33) 


3.3.2  Reduction  to  Quaternary  System 

Before  proceeding  to  develop  the  equations  for  six-  and  seven-component  sys¬ 
tems,  we  must  ensure  that  the  quinary  equation  will  reduce  to  that  of  a  quaternary 
system  under  the  proper  conditions.  This  is  accomplished  by  setting  z5  equal  to 
zero  and  K5  equal  to  one  (Warren,  1991). 

When  z5  becomes  zero,  so  must  x5  and  y5.  This  would  seem  to  leave  Af5 
undefined: 

lim  K<  =  —  »  undefined  (3.34) 

x3-*0  0 

We  can  remove  this  difficulty  by  the  application  of  1’HospitaTs  Rule.  The  expres¬ 
sion  becomes: 


lim  K5 

Xj— *0 
ys-*0 


dy5 

<&5_  _  l 
dx  5  1 

dx5 


(3.35) 


Therefore,  C5  =  Af5-l  =  l-  l=  0. 
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5 

To  avoid  division  by  zero,  multiply  Equation  (3.32)  by  ]~[  C,  : 

i=i 

a4  Ci  C2C3C4C5  + 

a3  £  (1  —  z  1)  C2C3C4C5  +  (i  —  z  2)  C 1  C3  C4  C5  +  (1  —  z  2)  C 1  C2  C4  C5  + 

(1  -  Z4)Ci  C2  C3  C5  +  (1  -  Z5)Ci  C2  C3  C4  j  + 

a2  (1  —  z  1  —  Z2)  C3  C4  C5  +  (1  —  Z}  —  Z3)  C2  C4  C5  +  (1  —  Zi  —  z4)  C2  C3  C5 

+  (1  -  2.1  —  Z5)  C2  C3  C4  +  (l  —  ^2  —  23)  C 1  C4C5  +  (1  ~ ■  Z2  —  Z4)  C 1  C3C5  + 

(1  -  Z2  —  Z5)  Cj  C3  C4  +  (1  — Z3  —  Z4)CiC2C5+(l  —  Z2  —  Z5)  Cj  C2C4  + 

(1  -  z4-z5)  Ci  C2C3j  + 

Ot  Z  1  C  i  (C  2  +  C3  +  C4+C5)  +  Z2C2(Ci  +  C3  +  C4  +  C5)  + 

Z3  C3  (C}  +  C2  +  C4  +  C5)  +  Z4C4(Ci  +  C2  +  C3  +  C§)  + 

Z  5  C  5  (C  1  +  C  2  +  C  3  +  C  4)  j  +  Z  i  Ci+Z2C2  +  Z3C3  +  Z4C4  +  Z5C5  =  0 

(3.36) 

Let  z5  and  C5  equal  zero: 
a3CiC2C3C4  + 

a2  J^(l  —  zi)C2C3C4+(l  — Z2)CiCiC4  +  (l  —  Z3)  Ci  C2C4  + 
(l-z4)CiC2C3]  + 

Ct  ^  Z  1  Ci  (C  2  +  C3  +  C4)  +  Z2C2  (C  i  +  C 3  +  C  4)  +  Z  3  C  3  (C  i  +  C 2  +  C  4)  + 

z4  C4  (C,  +  C2  +  C3)  +  Zi  Cj  +  Z2C2  +  z3C3  +  Z4C4  =  0  (3.37) 


This  is  identical  to  Equation  (3.21)  in  Warren’s  work,  which  is  the  expanded  form 
of  the  quaternary  equation.  We  can  now  safely  derive  the  expressions  for  six-  and 
seven-component  systems. 

3.3.3  Senary  System 


or  +  or 


'«  a -*.•>' 

+  a3 

f  f  (1  -  Zi  -  z^)* 

^  c 

.4=1  U1 

T  U 

.w  7=1+1  Q  C,  J 

or 


4  5  6  (1  -  Zi  ~  Zj  -  Zk) 

2*  2*  2j  c  r*  c 

i=l  j=i+ 1  *=i+  2  ci  '■'j  W 


a 


3  4  5  6  (1  -  Zj  ~Zj  -Zk-  Zf) 

I  Zj  Zj  Zj  Zj  C  f'  f' 

L*=I  7=1  + 1  4=1+2  1=1+3  W  W  c4  W 


6  Z; 


+  X  6 

n  c, 


=  0 


7*1 


(3.38) 


3.3.4  Septenary  System 


oc6  +  a5 


7  (1  -  Zj) 

h  Ci 


+  cr 


6  7  (1  -  2,  -  Zj) 

Zj  p  p 

i=1  7=1+1  <-i  W 


aJ 


5  6  7  (1  -  Z,  -  Z:  -  Zk) 

Z  Z  Z  -7T7 

i=1  7=«+l  4=i+2  W  L'4 


or 


4  5  6  7  (1  -  Z  -  Zj  -  Zk  -  Zt) 

I  Z  I  Z - r  FcT  — 

l_i=l  7=i+l  4=i+2  1=1+3  '-1  W  u*  W 


a 


Z  Z  Z  Z  X -(1  —  ”?L":r— z— — }- 

(.1=1  7=4+1  4=1+2  1=4+3  m=i+4 


C»  Cj  Ck  Ci  ~  Cm 


7  *i 

+  z  V—  = 0 


7 


(3.39) 
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3.4  Formulation  of  a  Generalized  Equation  for  the  Vapor  Fraction 


The  objective  function  can  be  recast  in  the  following  form: 

N  N 

X  Zj  C,  n  (1  +  a  Cj)  =  0,  where  C,,  Cj  =  constant  (3.40) 

i=  1  j*i 

Equation  (3.40)  is  in  the  form  of  the  generating  function  for  the  elementary  sym¬ 
metric  functions,  ar : 


HO +  '*«)=  I  tr  ar(Xi)  (3.41) 

i=l  r=0 

According  to  Macdonald  (1979),  a0(x,)  =  1  and  ar(xL)  =  0  for  all  r  >  n . 


We  can  now  express  the  objective  function  in  terms  of  the  r-th  elementary  sym¬ 
metric  function  in  C,  : 


N  N  N  N-l 

S  Zi  Ci  na+«  =  Z  zi  Ci  Z  °-r  ar(C  1 .  C« .  Cn)  =  0 

i=l  ]*i  i=l  r= 0 

where  C,  indicates  the  exclusion  of  the  i-th  term  from  the  operation. 
Since  ar  does  not  involve  i,  we  can  invert  the  order  of  the  summations: 


N- 1 

z  <*r 


r= 0 


-  £  z,-  C,  ar(C j,...,  C, CN)| 


=  0 


(3.42) 


(3.43) 


A  working  definition  of  the  elementary  symmetric  function  ar  could  be  "tak¬ 
ing  permutations  of  the  elements  of  a  set  r  terms  at  a  time."  For  example, 


fll(Cl,  Cj,..,,  Q)  =  (C1  +  C2+  •••  +■  Ctf) 


(3.44) 


35 


a2(C 


f1,  c2~\...,  cN-')  = 


CiC 


1'-'  2 


c,c 


lu3 


CiV-l  CjV 


(3.45) 


The  condition  C,  is  equivalent  to  the  j*i  condition  imposed  on  the  summa¬ 
tion  temis  in  the  earlier  versions  of  the  a-polynomial  and  herein  lies  the  computa¬ 
tional  awkwardness.  We  want  to  find  an  expression  that  allows  the  summation  to 
proceed  over  all  N  components,  which  is  an  operation  readily  represented  by  a 
DO-loop  in  computer  programming. 

To  eliminate  C,- ,  we  must  expand  the  symmetric  function.  In  Chapter  4,  we 
will  tackle  this  problem  after  a  discussion  of  symmetric  functions. 


Chapter  4 


DEVELOPMENT  OF  THE  GENERALIZED  EQUATION 


In  this  chapter,  we  shall  present  a  brief  introduction  to  the  theory  of  symmetric 
functions  to  show  why  they  provide  such  a  powerful  tool  to  express  permutations. 
Then  we  will  show  the  reasoning  used  in  the  search  for  a  recursive  expression  for  a 
in  terms  of  N ,  C,  and  zi .  Finally,  we  will  present  a  generalized  multicomponent 
equation  for  the  vapor  fraction,  a,  that  is  compact  and  readily  programmed  on  a 
computer. 

4.1  Introduction  to  Symmetric  Functions 

4.1.1  Notation  and  Definitions  of  Partitions 

Any  collection  of  v  non-negative  integers  (excluding  zero)  whose  sum  is  w  is 
called  a  v-partition  of  w.  The  individual  integers  are  referred  to  as  parts  of  the  par¬ 
tition  and  are  conventionally  written  in  descending  order  of  magnitude. 

David  et  al.  (1966)  state  that  if  there  are  X  distinct  parts,  say  P\,P2>—>P\ 
with  P\>P2>Pj>  •  •  •  >p\  £  1  and  if  p,  is  repeated  tt,  times,  with  1*1,  2,...,  X, 

then  the  partition  is  written  (pi'p*1  •  •  ■  Px1).  The  weight,  w,  of  the  partition  is 
written  as 

X 

W  =  £  PtKt 
i= 1 


(4.1) 
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and  the  number  of  parts,  or  length,  is 

v  =  2  *i  (4.2) 

i= 1 

Macdonald  (1979)  refers  to  jc,  as  the  multiplicity  of  i  in  the  partition.  For  example, 
the  partition  (42213)  has  weight  13,  6  parts  and  3  distinct  parts.  In  our  notation, 

Pi  =  4  and  Jtj  =  2;  p2  =  2  and  7^2  =  1 ;  p$  =  1  and  rc3  =  3 
4.1.2  Symmetric  Functions 


A  symmetric  function  is  one  in  which  the  individual  parts  can  be  interchanged 
without  altering  the  value  of  the  function,  such  as 

2*«  =  Xi+x2  +  x-}+  ■■■  +xn  (4.3) 

i=l 

The  number,  n,  of  the  quantities  x  does  not  affect  the  relationships  between  the 
various  forms  of  the  symmetric  functions,  but  does  appear  in  the  final  expressions. 
David  et  al.  (1966)  write 

2  Xi  =  0),  2  xi  ~  <r>  and  X  xixj  =  (")»  for  r  *  s  (4.4) 

i=l  1=1  i*j 

This  leads  directly  to  the  definitions  of  two  special  forms  of  symmetric  functions. 
MacMahon  (1920)  defines  the  unitary  or  a-functions  as 

ar  =  (  T  )  =  2  •  r  ~  1.  2,  •  •  •  (4.5) 

and  the  power  sums,  or  ^-functions,  as 

sr  =  (r)  =  2  xi  -  r  =  2>  '  '  ' 

i=l 


(4.6) 
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A  special  case  of  the  a-function  is  the  augmented  unitary  symmetric  function,  ur 
(David  et  al.,  1966): 

ur  =  [  lr  ]  =  r!(  V  )  =  r\ar  =  ,  (4.7) 

summed  over  all  ordered  sets  ix,  ...  ,  ir. 

4.1.3  Recursive  Expressions  for  Symmetric  Functions 

4.1.3.1  Interexpressibility  Tables 


Roe  (1918)  compiled  comprehensive  interexpressibility  tables  relating  the  vari¬ 
ous  classes  of  symmetric  functions  to  one  another.  These  consist  of  a  matrix  of 
coefficients  to  be  used  in  a  polynomial  which  might  yield,  for  example,  ur  =  /  (sr ). 
Of  interest  to  this  work  is  her  relationship  between  the  a-functions  (often  called  ele¬ 
mentary  functions)  and  the  5-functions;  it  is  presented  here  in  a  form  more  clearly 
expressed  by  David  et  al.  (1966): 


ar 


X  (-1  )('+»»> 

P?...P£ 


(4.8) 


David  et  al.  (1966)  also  used  this  equation  to  construct  interexpressibility 
tabl's  describing  polynomials  in  power-sum  series  (5)  for  a-functions  up  to  and 
including  weights  of  12.  F01  instance,  a  unitary  symmetric  function  of  weight  3 
would  be  represented  by  the  following  polynomial  from  their  Table  1.5.3: 


a3  =  -jr[(1)3  “  3(2)(1)  +  2(3) 


(4.9) 
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which,  in  terms  of  ^-functions,  is: 


a3  = 


3  S  2^j  +  2 s  3 


and,  when  written  as  power  sums,  becomes: 


(4.10) 


=  TT 


*  > 

n 

Z  *i 

. 

i=i 

3Z  xi2  Z  xi  +  2£  x? 

i=i  «=i  i=i 


(4.11) 


However,  neither  Equation  (4.8)  nor  Equation  (4.11)  is  conducive  to  solution  by 
computer  without  a  tremendous  table  look-up  effort. 


4.1.3.2  Determinant  Form 


Fortunately,  David  et  al.  (1966)  present  another  relationship  between  ar  and 
sr  in  determinant  form: 


ar  = 


1 

77“' 


i  . 

— -det 
r! 


*i 

1 

0 

0  • 

•  0 

s2 

*1 

2 

0  • 

•  0 

S3 

s2 

*1 

3  • 

•  0 

det  M 

Sa 

S  3 

^2 

^1  • 

•  0 

r! 

• 

. 

. 

• 

•  r-1 

Sr 

Sr-l 

sr-  2 

Sr -3  ' 

•  *i 

(4.12) 


This  provides  a  practical  method  of  calculating  ar  that  is  also  readily  programm¬ 
able  on  a  computer. 
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4.2  Search  for  a  Recursive  Expression  for  the  Vapor  Fraction 

Armed  with  a  working  knowledge  of  symmetric  functions,  we  can  manipulate 
the  a  equation  developed  in  Chapter  3: 


N-l  N 

X  ctr  *  X  zi  Cj  (C  !»•••*  C,-,...,  CN) 

r =0  1=1 


Y  =  o 


(3.43) 


to  eliminate  the  exclusion  term  C,  and  expand  the  symmetric  functions  into  a  more 
recognizable  form.  We  will  examine  the  results  for  several  values  of  r  and  use 
them  to  write  a  general  expression  for  a  as  a  function  of  N. 

4.2.1  Case  I;  r  =  N-2 

Equation  (3.43)  yields  the  following  coefficient  for  a: 


a<"-2>  J 


N 

X  zi  C,  aN-j(C  1  ,...,  C, ) 
1=1 


(4.13) 


We  can  expand  the  symmetric  function  a^_2  as  shown  in  Equation  (4.14).  Since 
the  exclusion  of  C,  from  the  product  on  the  RHS  gives  (N-\)  terms,  we  must  sub¬ 
tract  -J-  from  the  sum  to  yield  (N-2) : 
w 


<ty_2(C  i ,...,  Cj,...,  C^)  —  (Cv  -Cr  CN) 


1  1  c 
Cl  +  C2  + 


'N 


C, 


(4.14) 


To  eliminate  C, ,  we  can  write  the  product  on  the  RHS  of  Equation  (4.14)  as 

.  (Cy-Ctf) 

(Cv-Cr-CN)  =  — 77-^ 


(4.15) 


This  maneuver  will  allow  the  summation  to  proceed  over  all  N  components. 


After  substituting  Equation  (4.15)  into  Equation  (4.14),  we  have: 


&N— 2(C  Q CjV  ) 


(C \—CN) 

1 

,  1  .. 

1 

1  —  _ 

1 

Ci 

C,  1 

C2 

CN 

C, 

(4.16) 


Substituting  Equation  (4.16)  into  Equation  (4.13),  cancelling  C, ,  multiplying  by  z,  , 
and  then  summing  over  i  gives: 


0^-2)] 


(CrCN) 


'n  " 

I  z« 

«=i 


-L  +  J-  + 

Cl  c2 


1 

CN 


■if 

1=1 


(4.17) 


/v  n 

We  recall  that  £  z,  =  1  and  recognize  that  (Cf-Q)  =  ]""[  Q.  Noting  the  pres- 

i=l  *=1 


ence  of  an  elementary  symmetric  function  in 
Equation  (4.17)  as: 


-T+  ■■■  +  ± 

Ll  ln 


,  we  can  write 


a<*-2N 


N 

n  ck 

*= i 


i  i 
<V Cyv 


*  Zi 

-If 

1=1 


(4.18) 


4.2.2  Case  II:  r  *  N-3 


Equation  (3.43)  now  becomes: 


,(V-3)  J 


L  zi  C,  0^-3  (C Cl,...,  CN) 

i=l 


(4.19) 


We  can  expand  the  symmetric  function  %_3  as  shown  in  Equation  (4.20).  We 

.  1 
eliminate  C,  in  the  same  manner  as  in  Equation  (4.15)  and  remove  —  in  a  similar 

w 

fashion.  But  this  also  deletes  the  term  — which  is  necessary  to  cancel  the 
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C,-s 


corresponding  term  in  the  product.  Therefore,  we  must  compensate  by  adding  — 


%_3(Ci,...,  Ci,...,  Cfj)  - 


(C  i  —cN) 

1 

Ci 

CiC2 

1 

1 

p 

_I_  + 

1 

—  1 

...  +  J_ 

+  _l 

1 

*— » 

1 

'  C, 

.Cl 

c2 

c* 

J 

c,2. 

(4.20) 


Substituting  Equation  (4.20)  into  Equation  (4.19)  and  making  consolidations  similar 
to  the  previous  development  yields: 


a<*-3>  l 


N 

n  c* 

*=i 


i  i 

Cl  ’ cN 


-a  i 


1 


1 

CN 


N  Z:  N  2: 

Y  +  y  _ i_ 

*-•  r  ^  r  2 

,=l  ,=i  C, 


(4.21) 


4.2.3  Case  HI:  r  =  0 


We  have  saved  consideration  of  this  case  for  last  because  the  properties  of  a0 
are  not  readily  apparent.  It  would  seem  reasonable  to  interpret  a0(C],...,  C,,...,  CN) 
as  meaning  "taking  permutations  of  the  elements  of  a  set  zero  terms  at  a  time." 
However,  when  r  =  0,  ar  ->1  and  we  know  from  previous  developments  that  our 
a-polynomial  does  have  a  constant  term.  Therefore,  a0(C i,...,  C,,...,  CN)  must 
equal  one,  after  Macdonald  (1979).  So,  for  r  =  0,  Equation  (3.44)  becomes: 


N 

I  hCi 

»=1 


(4.24) 
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We  show  Equation  (4.18)  and  Equation  (4.21)  again  to  look  for  patterns  that  may 
assist  us  in  writing  the  expression  for  the  (N  -p  )-th  term: 


a 


CV-2) 


N 

n  ck 

k= 1 


Ol 


1  1 

Cl . cN 


N  Zi 

-if 

»=i 


(4.18) 


a' 


(N-  3) 


N 

n  ck 

k  = 1 


1  1 

C,  CN 


-a  i 


1  1 

Ci”"’  C* 


N  Zi  N  Zi 

i=1  i=l  W 


(4.21) 


4.2.4  The  General  Case:  r  =  (N-p),  p  =  1,2 . N 


After  substituting  for  r.  Equation  (3.43)  becomes: 


N 

£  a^N~p) 

p- 1 


N 

X  Zi  C|  flyv-p  (C]».”»  C;,...,  Cyv)  > 


=  0 


(4.23) 


By  continuing  the  expansion  of  this  equation  in  the  same  fashion  as  in  the  first  two 
cases,  we  note  a  descending  order  of  the  symmetric  function  and  an  ascending 
exponent  of  C,  with  each  additional  term.  This  leads  to  a  general  expression: 


N 
p- 1 


|  X  Zi  (C,-C*)[ap_i(C  j  1  Cyv1)  - 


Crlap_2(C Cy,-1)  +  C-2up_3(Cr1 . C*1)  -  ,...,  C£') 


±  .  C/71)  ±  Cy-^-1) 


=  0 


(4.24) 
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N 

Multiply  by  z,  ,  sum  i  from  1  to  N  and  recall  that  (C  j—  CN )  -  IT  '■ 

k= l 


N 

£  a{N~p) 
p= i 


'  N 

n  c* 
*=1 


w  z,  z,- 

flp-l  “  ap-  2  X  +  ap-  3  X  _  2  ~ 

»=1  i=l  c, 


+ 


JV 


Z; 


a  c/>-» 


=  0 


(4.25) 


N 

Since  n  Ck  does  not  involve  p,  we  can  move  this  term  outside  the  summa- 
*=i 

tion  sign  and  then  divide  it  out  as  a  factor  common  to  all  powers  of  a.  By  examin¬ 
ing  the  relationship  between  p,  the  subscripts  of  a  and  the  superscripts  of  Ci ,  we 
can  collapse  Equation  (4.25)  into  a  more  compact  form: 


N 


Z  a(*~p) 

p- 1 


i=l  '•'i 


=  0 


(4.26) 


where  ap_j  =  ap_j(  Cx  1 . C* 1  ) 


(4.27) 


a  0=  1 


(4.28) 


Ci=(Ki)T.F-  1 


(4.29) 
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Using  the  determinant  expression  for  the  elementary  symmetric  functions  that  was 
presented  in  4.1.3.2,  Equation  (4.27)  becomes: 

a,.;  »a,.>(Cf1 . <V  )  =  (4.30) 

The  matrix  M  has  dimensions  ip—j)  x  (p-j )  and  elements  given  by: 

* 

Sk-t+i  if  l  <  k 

[mkJ]='k  if  /  =  k+l  (4.31) 

0  if/>*+l 

n  r  i  y 

The  s  elements  are  given  by  sk  =  £  2,  (p-j)  (4.32) 

i=l  Ci 


Chapter  5 


VALIDATION  OF  THE  GENERALIZED  EQUATION 


The  first  test  of  validity  for  Equation  (4.26)  requires  that  it  be  equivalent  to  the 
form  of  the  objective  function  presented  in  Equation  (3.40).  Second,  it  must  gen¬ 
erate  the  same  coefficients  for  the  a  polynomial  that  were  produced  through  the 
expansion  of  the  objective  function  in  Equation  (3.25)  through  Equation  (3.32). 
Third,  the  equation  must  predict  the  proper  vapor  fraction  for  a  fluid  undergoing  an 
isothermal,  isobaric  flash  process. 

The  first  test  is  supplied  by  a  mathematical  proof  in  Appendix  A.  The  second 
test  can  be  accomplished  by  comparing  the  coefficients  produced  by  Equation  (4.26) 
with  those  of  Equation  (3.33).  Since  this  equation  has  already  been  shown  to 
reduce  to  that  for  a  quaternary  system  under  the  proper  constraints  on  z5  and  K5, 
then,  by  induction,  we  can  state  that  the  polynomial  produced  by  Equation  (4.26) 
will  do  the  same  and  therefore  should  be  valid  for  any  number  of  components. 

The  third  test  will  be  satisfied  by  comparing  the  equilibrium  ratios  generated 
by  Equation  (4.26)  with  experimental  values  determined  for  several  multicomponent 
hydrocarbon  fluids. 

5.1  The  Generalized  a  Equation  for  a  Quinary  System 

For  a  five-component  system,  Equation  (4.26)  becomes: 


\  f 

5  Z: 

(  iv+1  _  v  1 

■* 

11 

(  lr  ap_j  2,  j_x 

1  =  1  c  j 

> 

£  a(5_p) 

p=\ 


(5.1) 
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which  will  yield  a  quartic  polynomial  in  a: 

|i.4a4  +  p.3a3  +  \i2a2  +  ma  +  ^  =  0  (5.2) 

5.1.1  Coefficient  u4  (p  -  1) 


=  a0(Cj  !)  [  zj  +  z2  +  z3  +  z4  +  z5  ]  (5.3) 

We  have  already  said  that  a0(Q~!)  *s  defined  as  one  and  the  sum  of  the  mole  frac¬ 
tions  also  equals  one,  so  Equation  (5.3)  yields: 

1^4  =  1  (5.4) 


5.1.2  Coefficient  m  (p  ==  2) 


M-3  =  al (Ci  *)  S  2i  ~  a0 (Q  !)  Z  ~7T 
i= 1  i= 1  W 


^3 


1  —  z !  l-z2  1  -  r3  1  -  z4  1  -  z5 

Cj  c2  c3  c4  c5 


(5.5) 


(5.6) 


(5.7) 
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5.1.3  Coefficient  u?  (p  =  3) 


\h  =  a2(cr')  £  ,,  -  a,(Cr')  £  +  arfCf1)  £ 

i=l  i=l  Li  i=l  C,- 


(5.8) 


M-2  = 


1  +^  +  ^  +  ^  +  -J-  +  ^-  +  ^L-+  1 


CiC2  CjC3  c,c4  CjCs  c2c3  c2c4  C2c5  C3C,; 


1  1 
+  -  + 


^3^5  C4C5 


(D- 


11111 
C\  c2  c3  C4  C5 


21  _f2_  _^3_ 

Ci  c2  c3 


z4  z5 

+  c7+cT 


+  d) 


21  2  2  2  3  2  4  2  5 

+  - -  +  - -  +  - -  + 


C} 


C} 


cl 


cl 


cl 


(5.9) 


n  —  —  +  +  — — ~ —  +  +  +  — — — —  +  +  '”i' ~ + 

C1C2  CjC3  CtC4  C1C5  C2C3  C2C4  C2C5  C3C4 
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5.1.4  Coefficient  Hi  (p  =  4) 


Hi  =  a3(C-‘)£  2,  -  a2(Cf')i  jr  +  «,(Cf')i  -  a0<Cf'>i  7T  (5-12) 

i  =  1  «=1  W  1=1  Ci  i  =  1  C, 


Hi  = 


1  +  -1  +  — j — + —i_  +  — L_ + 


CiC2C3  C1C2C4  CxC2Cs  C1C3C4  C]C3C<, 


1  +  __L_  +  ^4^  +  _•  +  > 


C1C4C5  C2C3C4  C2C3C5  C2C4C5  C3C4C5 


(D- 


1  +_!_  +  ^+  ‘  +^V  +  ^L_  +  ^L-+  ' 


CjC2  C,C3  C,C4  CjCs  C2C3  C2C4  C2C5  C3C4 


C3C5  c4c5 


_2l_  J2_  ,  _f3_  _f4_  _*5_ 

Ci  c2  c3  c4  c5 


A_  +  JL  +  _L  +  _L  +  _L 

C !  C2  C3  C4  C5 


C2  /^2  /^2  /"•  2  /~  2 

1  L2  t-3  C4  C5 


(1) 


21  +  +  25 


c?  c,3  c,3  c? 


c<3 


(5.13) 


It  is  evident  that  the 


C,2C, 


terms  in  the  second  part  of  Equation  (5.13)  will  cancel 


those  in  the  third  part,  while  the  — -  terms  in  the  third  part  will  negate  the  entire 

Q 

fourth  part  of  the  equation.  The  first  and  second  parts  yield: 


1  -  Z,  -  Z2  -  Z3  1-Z1-Z2-Z4  1  -  Zi  -  Z2  -  Z5  I-Z1-Z3-Z4 

Hi  = - - - + - - - + - - - + 


C \C2C3 


C  ,C  2C  4 


C  iC  2C  5 


C1C3C4 
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1  -  Zl  -  z3  ~  z5  l-z1-z4-z5  1  -z2-z3  -z4  1  -z2-z3-z5 

+  - r  - + - —  r — r - + - - - r - + - r  .  _ - + 


C1C3C5 


C2C4C5 


C  2C3C  4 


1  -  Z2  -  Z4  -  Z^  1  -  z3  -  z4  -  z5 


c2cac5 


C3C4C5 


(5.14) 


5.1.5  Coefficient  Un  (p  =  5) 


Mo  -  «*<C,-1)  2  z,  -  »3(C,-')  i  -5-  +  a2(C-')  i-^7- 

«=1  i=l  L,- 


5 

z 

i=l 


<Ji(C,-‘)  2  7T  +  oo(cr’)  2  77 

«=1  I'i  i=l  W 


(5.15) 


The  analogous  cancellations  of  the  higher-order  -7-  terms  will  occur,  leaving  a 

w 


sum  of  five  terms  having  the  form 

1  ~  Z,  ~  Z,  ~  Zfc  ~  z< 

Q  Cj  Q 


(5.16) 


Since  the  mole  fractions  must  sum  to  one,  we  can  replace  the  numerator  of  Equa¬ 
tion  (5.16)  with  the  mole  fraction  of  the  remaining  component,  zm  ,  to  yield: 


Mo 


C  2C  3C  4C  5 


C  jC  3C  4C  5  C  \C  tC  &c 


2'-  4^  5 


z4  z5 

CiC2C3C5  C1C2C3CA 


(5.17) 


A  term-by-term  comparison  with  Equation  (3.32)  shows  that  the  generalized  a  poly¬ 
nomial  [Equation  (4.26)]  produces  identical  results. 
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5.2  Reproduction  of  Experimental  Vapor-Liquid  Equilibrium  Data 

5.2.1  Flash  Calculation  Package 

The  flash  calculation  package  used  in  this  work  incorporates  the  Af-value  equa¬ 
tion  of  Wilson  (1969)  and  the  modified  PREOS  proposed  by  Stryjek  and  Vera 
(1986a).  An  attempt  was  made  to  use  the  Af-value  prediction  of  Varotsis  (1989) 
but,  as  noted  in  Chapter  2,  it  was  developed  to  characterize  a  broad- spectrum 
petroleum  reservoir  condensate  or  crude  oil.  It  experiences  difficulty  handling  an 
arbitrary  hydrocarbon  mixture,  such  as  the  artificial  systems  for  which  equilibrium 
data  is  available  to  validate  this  work. 

The  volume  correction  of  Mathias  et  al.  (1989)  and  the  complementary  calcu¬ 
lation  of  Schick  and  Prausnitz  (1968)  for  mixture  pseudo-critical  volume  are  incor¬ 
porated  into  the  PRSV  EOS  but  since  it  is  only  required  to  generate  compressibility 
factors  and  fugacities,  the  modifications  have  no  noticeable  effect  on  the  computa¬ 
tions.  The  PRSV  EOS  shows  marked  improvement  over  the  PREOS  when  used  to 
duplicate  bubblepoint  and  dewpoint  studies  performed  by  Warren  (1991). 

The  binary  interaction  coefficients  used  in  the  PRSV  EOS  are  taken  from 
Knapp  et  al.  (1982)  and  Walas  (1985).  Physical  property  data  and  equation  param¬ 
eters  for  the  chemical  components  are  extracted  from  Stryjek  and  Vera  (1986b,c), 
Kumar  (1987)  and  Proust  and  Vera  (1989). 

The  computation  of  the  determinant  used  to  generate  the  elementary  symmetric 
functions  is  accomplished  with  a  modified  Gaussian  elimination  routine.  The  first 
elementary  symmetric  function,  alt  is  defined  by  a  [lxl]  matrix,  whose  deter¬ 
minant  constitutes  the  element  itself.  By  definidon,  a0  is  set  equal  to  one. 
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The  polynomial  is  evaluated  at  the  bubblepoint  line  (a  =  0)  and  an  interval¬ 
halving  technique  is  used  to  march  across  the  two-phase  region  until  the  value  of 
the  polynomial  changes  sign,  indicating  the  vicinity  of  the  root  Then  a  Newton- 
Raphson  iterative  search  is  conducted  to  converge  to  the  exact  value  of  a. 

5.2.2  Binary  System 

The  fugacity-based  flash  algorithm  is  used  to  replicate  the  equilibrium  ratios 
determined  by  Bloomer  et  al.  (1953)  for  a  methane-ethane  system  at  a  temperature 
of  -60  °F  over  a  pressure  range  of  100-900  psia.  Comparisons  of  calculated  and 
empirical  values  of  ATCH4  and  appear  in  Figure  5.1  and  5.2,  respectively.  The 
results  lie  within  the  margin  of  error  attributable  to  the  PRSV  EOS. 

5.2.3  Septenary  System 

Standing  (1977)  provides  a  sample  flash  calculation  for  a  seven-component 
hydrocarbon  system  reported  by  Dodson  and  Standing  (1941),  complete  with  values 
for  experimental  AT,  and  the  vapor  fraction.  This  sort  of  data  allows  the  calculation 
of  a  solely  on  the  basis  of  computing  the  coefficients  of  the  a-polynomial  and 
determining  the  applicable  root,  with  no  recourse  to  the  equation  of  state.  Once  the 
interval-halving  search  provides  an  initial  estimate  of  the  root,  the  Newton-Raphson 
technique  converges  in  three  iterations  to  a  value  of  a  identical  to  that  calculated  by 
Standing. 


Equilibri 
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5.2.4  Predicting  Roots  with  the  Fourier-Budan  Theorem 

A  useful  theorem  for  predicting  the  number  of  roots  of  a  polynomial  that  can 
occur  on  a  particular  interval  is  that  of  Fourier  and  Budan  (Barbeau,  1989).  Sup¬ 
pose  p(t)  is  a  polynomial  over  the  field  of  real  numbers,  R,  and  that  u  and  v  are  real 
numbers  with  u  <  v  and  p(u)p(v )  *  0.  The  number  of  zeros  between  u  and  v 
cannot  be  greater  than  A  -  B ,  where  A  is  the  number  of  changes  of  sign  in  the 
sequence  {  p(u),  p\u),  p"(u),  ...,  p^n\u)  )  and  B  is  the  number  of  changes  of 
sign  in  the  sequence  {  p(v),  p'(v),  p"(v ),...,  p(n)(v)  }■  If  this  number  differs 
from  A  -  B ,  it  must  do  so  by  an  even  amount 

An  interesting  aspect  of  the  polynomial  expression  for  the  vapor  fraction  is  its 
capability  to  mathematically  confirm  the  existence  of  a  unique  value  within  the 
two-phase  region  for  a  particular  set  of  feed  conditions.  This  is  equivalent  to  stat¬ 
ing  that  the  polynomial  has  only  one  zero  on  the  interval  0  <  a  <  1.  From  the  phy¬ 
sics  of  the  problem,  we  know  this  to  be  true  but,  by  the  use  of  the  Fourier-Budan 
theorem,  we  can  also  prove  it  rigorously. 

Let  us  test  this  theorem  on  the  septenary  system  of  Standing  (1977)  utilized  in 
5.2.3;  this  is  represented  by  a  sixth-order  polynomial; 

ma6  +  n5a5  +  ii4a4  +  n3a3  +  \i2a2  +  m<x  +  Po  =  0  (5. 18) 

where 


m>  = 

-9.58519 

H3  = 

87.24949 

hi  = 

65.90501 

H4  * 

-21.71701 

H2  = 

-120.72959 

H5  = 

-1.76522 

1.00000 
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We  can  differentiate  Equation  (5.18)  six  times  and  form  the  derivative  sequences 
for  u  =  0  and  v  =  1.  The  sign  changes  are  summarized  in  Table  5.1. 
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Since  A  -  B  =1,  there  exists  only  a  single  root  of  the  polynomial  on  the 
interval  [0,1];  an  examination  of  the  graph  of  the  function  (Figure  5.3)  confirms  this 
fact.  Therefore,  we  can  use  the  interval-halving  and  Newton-Raphson  solution  pro¬ 
cedures  outlined  at  the  beginning  of  this  chapter  with  confidence  that  they  will 
obtain  a  unique,  realistic  value  of  the  vapor  fraction. 

5.2.5  Decenary  System 

Gregory  et  al.  (1971)  performed  equilibrium  measurements  on  a  lean  natural 
gas  at  cryogenic  conditions.  It  is  reported  as  a  ten-component  system  with  the  feed 
composition  shown  in  Table  5.2.  This  is  a  very  "sparse"  ten-component  gas,  with 
six  components  present  in  extremely  dilute  concentrations.  The  ^-values  for  the 
last  four  constituents  were  zero  for  eleven  of  the  sixteen  operating  conditions  tested 
in  this  work,  denoted  in  Table  5.3  by  the  run  number  assigned  by  the  investigators. 
The  remaining  twelve  sets  of  published  data  duplicate  conditions  in  one  of  the 
tested  runs  or  are  incomplete  due  to  apparatus  failure.  The  use  of  the  Fourier- 
Budan  theorem  provides  warning  that  perhaps  this  gas  would  be  better  represented 
by  an  equivalent  "lumped"  system. 

Recall  that  the  number  of  roots  predicted  by  the  Fourier-Budan  theorem  is  the 
maximum  possible  and  may  differ  from  the  true  value  by  only  an  even  integer. 
This  is  demonstrated  in  Table  5.4,  where  both  the  predicted  and  actual  number  of 
roots  for  each  run  are  tabulated.  The  Newton-Raphson  technique  converges  to  the 
experimental  value  for  ten  of  the  sixteen  runs.  Three  other  data  points  follow  the 
proper  trend,  while  no  root  is  found  on  the  interval  [0,1]  for  three  other  conditions 
(Figure  5.4). 


Table  5.2  -  Feed  Composition: 
10'Component  Natural  Gas 
(Gregory  et  al.,  1971) 

Component 

2, 

Component 

Nitrogen 

0.00600 

n-Butane 

0.00070 

Methane 

0.95970 

/-Pentane 

0.00030 

Ethane 

0.03000 

n-Pentane 

0.00010 

Propane 

0.00390 

3-Methylpentar.e 

0.00025 

/-Butane 

0.00070 

2-Methylhexane 

0.00015 

Table  5.3  -  Experimental  Flash  Conditions: 
10-Component  Natural  Gas 
(Gregory  et  al.,  1971) 


Run 

Pressure 

(psia) 

Temperature 

(°F) 

Run 

Pressure 

(psia) 

Temperature 

(°F) 

1 

300.0 

-156.3 

14 

100.0 

-200.0 

3 

100.0 

-206.0 

15 

500.6 

-127.0 

B 

700.0 

-103.0 

18 

23.0 

-252.0 

B 

500.0 

-125.0 

20 

497.0 

-129.0 

8 

498.5 

-120.0 

21 

23.5 

-251.5 

9 

695.0 

-105.0 

25 

700.0 

-107.0 

10 

100.0 

-203.3 

26 

298.0 

-157.5 

12 

100.0 

-195.0 

28 

500.0 

-130.0 

Table  5.4  -  Results  of  the  a- Polynomial  and  Fourier-Budan  Theorem: 
10-Component  Natural  Gas 
(Gregory  et  al.,  1971) 
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An  examination  of  graphs  of  the  polynomial’s  behavior  over  a  range  of  a  for 
Runs  3  and  4  (Figures  5.5  and  5.6)  confirms  the  algorithm’s  prediction  that  no  roots 
exist  within  the  phase  envelope.  The  case  of  Run  9  is  not  so  obvious.  Its  graph 
(Figure  5.7)  shows  that  the  function  exists  entirely  above  the  abscissa;  hence  no 
root  is  possible.  However,  if  the  resolution  of  the  graph  is  increased  to  examine  the 
region  very  near  the  axis,  two  local  minima  are  revealed  (Figure  5.8).  One  of  these 
corresponds  to  the  experimental  value  of  a  determined  for  this  run.  The  polyno¬ 
mial  is  attempting  to  represent  the  system’s  behavior  but  is  not  completely  success¬ 
ful  because  the  low  concentration  of  certain  components  effectively  prevents  the  gas 
from  acting  like  a  decenary  system. 

It  is  instructive  to  compare  the  form  of  the  a-polynomial  with  that  of  the 
Rachford-Rice  objective  function  which  is  superimposed  on  Figure  5.7.  The  same 
high-resolution  scan  of  the  graph  of  the  latter  equation  depicts  no  equivalent  max¬ 
ima  which  might  identify  the  vapor  fraction  in  the  manner  of  the  polynomial. 

5.2.6  Lumping  a  Decenary  System  into  a  Quaternary  System 

The  a-polynomial  successfully  converges  to  the  proper  answer  for  a  majority 
of  the  runs;  however,  it  also  yields  multiple  roots  where  the  physics  of  the  problem 
allows  only  one.  This  suggests  that  the  system  is  not  being  properly  modeled.  The 
categorization  of  the  fluid  as  a  ten-component  natural  gas  is  overly  generous  in  light 
of  the  fact  that  six  of  its  chemical  constituents  are  present  in  mole  fractions  meas¬ 
ured  in  the  ten-thousandths.  It  was  decided  to  represent  this  sparse  gas  as  a  four- 
component  lumped  system,  consisting  of  methane,  ethane,  nitrogen  and  propane. 

The  mole  fractions  of  this  new  fluid  are  normalized  and  the  resulting  cubic 
polynomial  in  a  is  solved.  The  Fourier-Budan  theorem  predicts  a  maximum  of  one 
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root  in  the  two-phase  region,  to  which  all  sixteen  runs  converge.  The  numbers 
tabulated  in  Table  5.5  and  displayed  graphically  in  Figure  5.9  attest  to  the  validity 
of  this  lumping  scheme.  An  attempt  was  made  to  eliminate  the  next  leanest 
component— propane-- from  the  mixture  and  model  the  system  as  a  ternary,  but  this 
resulted  in  spurious  roots  for  all  data  runs  and  was  hence  rejected  as  unrealistic. 


Table  5.5  -  Results  of  the  a-Polynomial  and  Fourier-Budan  Theorem: 
"Lumped"  4-Component  Natural  Gas 
(Gregory  et  ah,  1971) 


Run 

Root  Limit 
(Actual) 

Newton- 

Raphson 

Iterations 

Roots  on  [0,1] 

Initial  Guess 

Calculated 

Experimental 

1 

1 

3 

0.635 

0.631 

0.603 

3 

1 

3 

0.035 

0.038 

0.155 

D 

1 

3 

0.995 

0.998 

0.911 

D 

1 

3 

0.815 

0.819 

0.761 

8 

1 

3 

0.935 

0.932 

0.908 

9 

1 

3 

0.925 

0.920 

0.795 

10 

1 

3 

0.715 

0.711 

0.687 

12 

1 

3 

0.895 

0.899 

0.890 

14 

1 

3 

0.845 

0.842 

0.830 

15 

1 

3 

0.765 

0.770 

0.747 

18 

1 

3 

0.035 

0.037 

0.109 

20 

1 

3 

0.635 

0.631 

0.591 

21 

1 

3 

0.015 

0.012 

0.078 

25 

1 

3 

0.775 

0.779 

0.548 

26 

1 

3 

0.475 

0.473 

0.430 

28 

1 

3 

0.535 

0.536 

0.486 

Chapter  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Conclusions 

1.  The  Rachford-Rice  objective  function  can  be  represented  as  a  polynomi  '  in  a, 
the  system  vapor  fraction.  Its  coefficients  involve  elementary  symmetric  func¬ 
tions,  which  can  be  expressed  in  terms  of  a  determinant  whose  elements  are 
functions  of  equilibrium  ratios  and  feed  composition. 

2.  The  polynomial  has  been  shown  to  be  well-behaved  within  the  two-phase 
vapor-liquid  region  if  the  system  is  properly  defined  in  terms  of  the  number  of 
its  components.  The  vapor  fraction  root  on  the  interval  [0,1]  can  be  quickly 
determined  using  an  ordinary  interval-halving  technique  to  provide  an  initial 
estimate  to  the  Newton-Raphson  iterative  method. 

3.  The  regular  behavior  of  the  polynomial  lends  itself  to  use  as  a  descriptive  tool 
for  the  conditions  of  the  system  within  the  phase  envelope.  The  Rachford- 
Rice  objective  function  is  not  capable  of  this  task  as  evidenced  by  Figure  5.7; 
its  unpredictable,  singular  nature  offers  no  clue  to  the  reason  a  root  was  not 
found  on  the  interval  [0,1]  for  this  case.  As  discussed  earlier,  a  close  exami¬ 
nation  of  the  curve  of  the  polynomial  revealed  a  local  minimum  at  the  experi¬ 
mental  value  of  a.  This  became  a  realistic  root  (a  <  1)  once  the  system  was 
lumped  into  four  components. 
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4.  The  theory  of  polynomials  is  well-developed  and  their  behavior  and  zeros  can 
be  predicted  with  good  confidence.  By  the  use  of  the  Fourier-Budan  theorem, 
it  can  be  shown  mathematically  that  only  one  real  root  for  the  a-polynomial 
can  exist  on  the  interval  [0,1]  for  a  well-defined  system.  This  eliminates  the 
need  to  solve  for  all  the  roots  of  an  N-\h  order  polynomial. 

5.  The  Fourier-Budan  theorem  can  be  used  as  a  tool  for  investigating  various 
lumping  schemes  whereby  multicomponent  fluids  are  condensed  to  equivalent 
systems  composed  of  fewer  components.  The  phase  behavior  of  sparse  fluids 
having  dilute  concentrations  of  several  constituents  does  not  seem  to  be  well- 
described  by  the  polynomial  of  degree  appropriate  to  the  number  of  com¬ 
ponents.  In  this  case,  the  polynomial  yields  no  roots  or  at  least  two  roots 
inside  the  phase  envelope  for  certain  temperature  and  pressure  conditions.  It 
appears  that  a  lumping  scheme  can  be  tuned  by  generating  pseudocomponents 
to  give  successive  polynomials  of  lower  degree  until  only  one  root  is  deter¬ 
mined  on  the  interval  0  <  a  <,  1. 

6.2  Recommendations 

1.  Further  study  should  focus  on  coupling  the  polynomial  algorithm  to  an  equa¬ 
tion  of  state  and  extending  this  work  to  flash  calculations  involving  more  than 
two  phases. 

2.  Timing  studies  could  be  conducted  to  determine  the  exact  savings  in  CPU  time 
realized  by  the  use  of  the  polynomial  instead  of  the  Rachford-Rice  objective 
function. 
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3.  Peng  et  al.  (1975)  estimate  that  75%  of  the  total  computing  time  in  composi¬ 
tional  reservoir  simulation  may  be  related  to  the  phase-behavior  part  of  the 
program.  The  savings  in  computational  workload  realized  by  the  use  of  the 
generalized  equation  developed  in  this  work  might  be  applied  to  the  employ¬ 
ment  of  an  EOS  better  able  to  predict  fluid  thermodynamic  properties.  The 
highly  nonlinear  nature  of  the  equations  proposed  by  Benedict,  Webb  and 
Rubin  (1940,  1942,  1951)  or  Lee  and  Kessler  (1975)  require  iterative  solutions 
but  they  yield  much  more  accurate  representations  of  fluid  behavior,  espe¬ 
cially  of  nonhydrocarbon  systems. 

4.  Since  the  coefficients  of  the  generalized  polynomial  depend  only  on  the  feed 
composition  and  equilibrium  ratios,  research  should  continue  to  develop  highly 
accurate  K-value  prediction  methods  (e.g.,  on  the  basis  of  convergence  pres¬ 
sure).  If  this  can  be  done  with  sufficient  accuracy,  the  fugacity-convergence 
approach  and  its  inherent  dependence  on  an  equation  of  state  can  be  sup¬ 
planted  for  flash  calculations  where  nothing  more  than  the  phase  split  and 
compositions  are  required.  The  polynomial  algorithm  can  be  solved  on  a  pro¬ 
grammable  scientific  calculator  and  would  provide  the  engineer  with  a  valuable 
predictive  tool  in  situations  where  he  or  she  has  no  access  to  a  computer  capa¬ 
ble  of  running  an  EOS-based  flash  routine. 
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Appendix  A 

MATHEMATICAL  PROOF  OF  THE  GENERALIZED  EQUATION 
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KNOWN: 


X  *<  Ci  n  (1  +  «  Cj)  =  £  a^>  X  z,  C,  J  a^C, . C,,...,  CN)\  (A-l) 

»=i  y*  p=i  i=i  l  J 

POSTULATE: 


*  N  N  ,K1  s  *  if  N  Z:  1  1 

I  ^  n  a  +  «  Cj)  =  £  a!»-»  j  n  c*  £  <-iy*‘  a,.;  £  — ^  U-2) 

i= 1  j*i  p  =  1  [  4=1  _/=l  _  i=l  L,  J 


=  £  a<"T»  £  z,C,  |  n  Ct  £  (-1V*1  HA-3) 

p=i  i=i  [4=1  y=i  I  C,  J 


Prove  that  the  coefficients  of  a  in  Equation  (A-l)  and  Equation  (A-3)  are 


equivalent: 


«W-,(C  1 C, c„)l  =  In  c4  £  (-iy«  Seu 

J  14=1  j= 1  Q' 


PROOF: 


(A-4) 


We  can  express  the  a-function  as: 

N 

n  ct 

ON-piCi . Ci,...,CN)  =  ~ —  ap_,(Cr‘ . Cf' . Cf1)  (A-5) 


A/ 

nck 

where  — — —  represents  (N- 1)  terms:  N-p  =  (N- 1)  -  ip- 1) 

Lj 

Eliminate  the  C,-1  term  in  the  RHS  of  Equation  (A-5)  by  rewriting  the  a-function 


flp-i (Cf1 . C,1 . C*1)  =  ap.^Ci1 . <V)  -  Cf1  flp_2(Cr1,...,CArI1)(A-6) 


-1  /°.-i'i 


-1  r  -1 
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ap-l(P\ 1  ••••>  C't  1>—' »  CN*)  -  flp_2(Cf1 C/71)  -  C,  1  ap_3(C i 1 C/yii  )(A-7) 
flp-3(Cj 1 C,-  C/^1)  =  ap_3(Ci 1 C/71)  -  C,-  1  flp_4(C j 1 C/yli  )(A-8) 


^-(p-DCCf1 Q  1,..„  C/71)  -  flp-<p-i)(Cf 1 C/71)  - 

Cr^floCCf1 . C/7li)  (A-9) 

Recall  that  a0  =  1  and  then  substitute  Equations  (A-7),  (A-8),...,  (A-9)  into  Equation 
(A-6): 

flp.jCCj1,...,  C,  1 . O  =  flp-i(C ! 1 C#-1)  —  C,-  1  flp_2(C j 1 C/y1) 

+  C,  2  a/,_3(Cf1 C/71)  -  C,-3  flp-^CCf1 C/71)  + 

•  •  •  ±  C,~<^_2>  Qp-{p-i){C f1 C/71)  ±  cr^-x)  (A- 10) 

After  writing  the  recursive  form  for  the  RHS  of  Equation  (A- 10),  the  equation 
becomes: 


a,-i(C|-' . =  £  H r'C.-O-Ufl,., 

/=> 

Substitute  Equation  (A- 11)  into  Equation  (A-5): 

N 

nc, ,  -i 

aN-p(C, C, •’•••»  Cw)  =  £  (-iy+1 

L/=1 

Combine  C,  terms: 

. c,,...,  cN)  =  n  c*  i  (-iy+1 

*=i  L/=i  c,  y  J 


(A- 11) 


(A- 12) 


(A- 13) 


Equation  (A— 13)  =  Equation  (A-4)  Q.E.D. 


Appendix  B 

ALGORITHM  FLOWCHART 


Appendix  C 
COMPUTER  CODE 
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***************************************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 


♦ 

* 


7  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 


Advisor 


M.S.  thesis 

Dr.  Michael  A.  Adewumi 


******************************************************************* 

*  Program  ALPHATEST  (FORTRAN  77) 
******************************************************************* 

* 

*  This  program  calculates  values  of  the  vapor  fraction, 

*  given  equilibrium  ratios,  Ki,  and  feed  mole  fractions, 

*  zi.  It  can  be  used  to  reproduce  experimental  results  of 

*  equilibrium  flashes. 

* 

*  ALPHATEST  calls  ALPHACOEFF.  BUD  AN,  ALPHAPLOT,  and 

*  ALPHAROOT 

*  ALPHACOEFF  calls  subroutine  SYMFUNCTION 

*  SYMFUNCTION  calls  subroutine  DETERM  and  function  FACTOR 


VARIABLES: 


alpha  =  calculated  system  vapor  fraction 
beta  =  experimental  system  liquid  fraction 
coefficient  =  coefficient  of  alpha  polynomial 
Ki  =  equilibrium  ratio  for  component  i 
molefrac  =  feed  mole  fraction  of  component  i 
Ncomp  =  number  of  components  in  feed 
Npress  =  number  of  data  sets  to  be  evaluated 
Pi  =  system  pressure,  psia 
Ti  =  system  temperature,  F 
xalpha  =  experimental  system  vapor  fraction 


It  is  formatted  to  input  zi,  temperature,  pressure,  liquid 
mole  fraction,  and  Ki 


********************************************************************* 


IMPLICIT  REAL*8(a-h,o-z) 

REAL*8  Ki(500,100)4nolefrac(0:100) 
PARAMETER(Npress=  1 6  ,Ncomp=  1 0) 

DIMENSION  alpha(SOO),  beta(500),  coefficient(0:100). 
@  Pi(500),  Ti(500),  tanay(2),  xalpha(500) 

Data  Input 


The  number  of  components  (Ncomp)  and  the  number  of  data  sets 
to  be  run  (Npress)  arc  specified  as  PARAMETERS' 


Open  and  Rewind  Input  and  Output  Files 

OPEN(unit=  1  ,file='indata',status=  'old”) 
OPEN(unit=7,file='table',status=  'unknown') 
OPEN(unit=8,file='plot',status='unknown') 

REWIND(unit=l) 

REWIND(unit=7) 

REWIND(unit=8) 

read(l/)  (molefrac(i),  i  =  1,  Ncomp) 
do  1000  j  =  1,  Npress 

read(l ,*)  Pig),  Tig),  betag) 
read(l,*)  (Ki(j,i),  i  =  1,  Ncomp) 
xalphag)  =  l.dO  -  betag) 

1000  continue 

Choose  between  single  or  multiple  runs 

write(6,*)  'Evaluate  one  data  set?  enter  V 
write(6,*)  'Evaluate  all  data  sets?  enter  2' 
read(5,*)  numsets 

if(numsets  .EQ.  1)  then 

write(6,*)  'Enter  number  of  data  set  for  this  run' 
read(5,*)  j 
go  to  2100 
end  if 

do  2000  j  =  1,  Npress 

2100  write(7,*)  '  ' 
write(7,*)  '  ' 

write(7,*)  '  RUN  'j 

write(6,*)  'J  =  ',  j 
write(7,2500)  Pig),Tig),betag) 

2500  formatCTressure  =  ',f6.1,'  psia  Temperature  =  ',f6.1,'  F 
(©Liquid  Mole  Fraction  =  ',f6.4) 


Call  subroutines 

Calculate  coefficients  of  polynomial 

call  ALPHACOEFFCNcompJ^pressjjnolefrac.Ki .coefficient) 


Predia  the  number  of  roots  on  [0,1]  by  Fourier-Budan  theorem 


call  BUDAN(j,Ncomp,coefficient,numroot) 
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Solve  for  the  roots  by  Newton-Raphson  method 

call  ALPHAROOT(j,Ncomp,coefficient,xalpha,numroot,alpha) 

Generate  various  plots  (EDIT  the  file  to  remove  comments  for  specific 
options) 


call  ALPHAPLOT (Ncom pjjnolefrac, alpha, coefficient.Ki) 

*  ALPHAROOT  has  internal  output  section  to  compile  a  table 

*  listing  statistics  on  the  determination  of  alpha 

2000  continue 

********************************************************************* 

*  Produce  this  format  to  plot  data  points  as  dots: 

*  (PLOTFAT=20) 

*  2 

*  x(l)  y(l) 

*  x(l)  y(l) 

*  2 

*  x(2)  y(2) 

*  x(2)  y(2) 

*  etc. 

********************************************************************* 


do  3000  j  =  1,  Npress 

write(8,3500)  alpha(j).xalpha(j),alpha(j),xalpha(j) 

3500  format('2  V.e  1 6.9, 10x,e  1 6.9  J,e  1 6.9,1  Ox  ,e  1 6.9) 

3000  continue 

CLOSE(unit=l) 

CLOSE(unit=7) 

CLOSE(unit=8) 

stop 

aid 

******************************************************************** 
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*  University  Park,  Pennsylvania 

* 

*  M.S.  thesis 

* 

*  SUBROUTINE  ALPHACOEFF 

* 

*  This  subroutine  calculates  the  coefficient  for  each  term  in 

*  the  general  polynomial  for  the  vapor  fraction,  alpha: 

* 

*  P(alpha)  =  cO  +  cl*alpha  +  c2*alpha**2  +  ...  + 

* 

*  cOsicomp-l)*alpha*,*,(Ncomp-l) 

*  Equation  4.29  in  the  thesis. 

* 

********************************************************************* 


SUBROUTINE  ALPHACOEFF(NcompJ4pressJjjnolefrac,Ki, coefficient) 

IMPLICIT  REAL*  8(a-h  ,o-z) 

REAL* 8  Ki(500,100),  molefrac(0:100) 

INTEGER  p 

DIMENSION  coefficient(0: 1 00),  c(100) 

OPEN(unit=  14,file=  'coefT,status=  'unknown ") 

OPEN(unit=  1 5  ,file= 'coeff. plot ',status= 'unknown  0 

if(Ncomp  .LT.  2)  then 

write(6,*)'You  cannot  flash  this  system' 
stop 
end  if 


* 

*  Calculate  Ci  =  Ki  -  1 

* 


do  0500  k  =  1,  Ncomp 
c(k)  =  Ki(jjjc)  -  l.dOO 
0500  continue 


*  p— loop  increments  the  power  of  alpha 

* 


C  write(15,*)Ncomp 

do  1000  p  =  1,  Ncomp 
temporary  =  O.dOO 

do  2000  j  =  1,  p 

*  Zero-order  elementary  symmetric  function,  a0[l/G],  defined  as  1 

if(p— j  .EQ.  0)  then 
apj  =  l.dOO 
go  to  2500 


n  n 


end  if 


Call  subroutine  to  calculate  the  elementary  symmetric 
function,  apj 

call  SYMFUNCTION (Ncomp  j,p,c,apj) 

ratio  =  O.dOO 

do  3000  i  =  1.  Ncomp 

ratio  =  ratio  +  molefrac(i)/(c(i)**(j-l)) 
continue 

temporary  =  temporary  +  ((-l.dO)**(j+l))*apj*ratio 
continue 

coefficient(Ncomp-p)  =  temporary 

write(14,*)'Coefficient(',Ncomi>-p,0  =  ',coefficient(Ncomp-p> 
write(15,*)Ncomp-p,coefficient(Ncomp-p) 

1000  continue 

return 

end 

*********************************************************** ******** 
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BRETT  D.  WEIGLE 
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M.S.  thesis 

SUBROUTINE  SYMFUNCTION 
This  subroutine  calculates  the  elementary  symmetric  function 

a(p— j){  1/Ci} 

*************************************** **************************** 


* 

* 

2500 

3000 

2000 


SUBROUTINE  SYMFUNCTION  (Ncomp  J,p,c,  apj) 

IMPLICIT  REAL*8(a-h,o-z) 

REAL*8  mmatrix(100,100) 

INTEGER  factor, p 


DIMENSION  c(100),  s(100) 


Compute  the  power-sum  series:  s  =  sigma[  (1/Ci)**lambda  ] 


n  =  0  -  j 

do  1000  lambda  =  1,  n 
sum  =  O.dOO 

do  2000  i  =  1,  Ncomp 

sum  =  sum  +  (l.d0/c(i))**lambda 
2000  continue 

sflambda)  =  sum 
1000  continue 

Build  the  matrix  MMATRIX 

do  3000  k  =  1.  n 

do  4000  1  =  1,  n 

if(l  .LE.  k)  mmatrix(kj)  =  s(k-l+l) 
ifQ  .EQ.  k+1)  mmatrixOc.l)  =  DFLOAT(k) 
if(l  .GT.  k+1)  mmatrix(k,I)  =  O.dOO 
4000  continue 

3000  continue 

Since  al{l/Ci}  forms  a  [1x1]  matrix,  its  determinant  is  the 
element  itself 

if(p-j  .EQ.  1)  then 
det  =  mmatrix(l.l) 
go  to  5000 
end  if 

Compute  the  determinant  of  MMATRIX 


call  DETERM(mmatrix,n,det) 


Compute  the  elementary  symmetric  function 


5000  apj  =  det/factor(n) 

return 

end 


Function  to  compute  the  factorial 


FUNCTION  factor(n) 

INTEGER  factor,!  ,n 

factor  =  1 
if(n  .GT.  0)  then 
do  6000  i  =  2jt 
factor  =  factor*i 
6000  continue 
end  if 
end 

******************************************************************** 
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SUBROUTINE  DETERM 

This  program  calculates  the  determinant  of  an  NxN  matrix. 

First,  partial  pivoting  is  performed  on  a  nonsingular  matrix  by 
Gaussian  elimination.  This  produces  a  triangular  matrix  whose 
determinant  can  be  calculated  by  computing  the  product  of  all 
the  diagonal  entries. 

The  augmented  matrix  does  not  contain  the  normal  last  column 
which  represents  the  right-hand  side  of  a  system  of  linear 
equations;  AUG  is  the  same  as  the  original  matrix. 

VARIABLES: 

N  =  dimension  of  matrix 
AUG  -  augmented  matrix 
U.K  =  indices 

MULT  =  multiplier  used  to  eliminate  an  unknown 
PIVOT  =  used  to  find  nonzero  diagonal  entry 

******************************************************************* 


SUBROUTINE  DETERM(augji.det) 

IMPLICIT  REAL*  8(a-h,o-z) 
REAL*  8  mult 
INTEGER  pivot 
DIMENSION  aug(100,100) 


* 

* 


Gaussian  elimination 


do  7000  i  =  1,  n 
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Locate  nonzero  entry 


if(aug(i,i)  -EQ.  0)  then 
pivot  =  0 
j  =  i  +  1 

3000  if((pivot  .EQ.  0)  .AND.  (j  .LE.  n))  then 

if(aug(j,i)  .NE.  0)  pivot  =  j 
j  =  j  +  1 
go  to  3000 
end  if 

if(pivot  .EQ.  0)  then 

print  *, 'Matrix  is  singular' 
stop 
else 


Interchange  rows  I  and  PIVOT 


do  4000  j  =  i,  n 
temp  =  aug(ij) 
aug(i,j)  =  aug(pivotj) 
aug(pivotj)  =  temp 
4000  continue 

end  if 

end  if 


Eliminate  I— th  unknown  from  equations  1+1,  ...,  N 


do  6000  j  =  i+1,  n 

mult  =  -aug(j4)  /  aug(i.i) 

do  5000  k  =  i,  n 

aug(j,k)  =  aug(j,k)  +  mult  *  aug(ijc) 

5000  continue 

6000  continue 

7000  continue 

Calculate  the  determinant  of  matrix  AUG  by  computing  the 
product  of  the  diagonal  elements 

prod  =  l.dO 
do  8000  i  =  1,  n 
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do  9000  j  =  1,  n 

if(i  .EQ.  j)  prod  =  prod  *  aug(i,j) 

9000  continue 

8000  continue 

det  =  prod 

return 

end 

******************************************************************** 

4  Dec  91 

BRETT  D.  WEIGLE 

Petroleum  and  Natural  Gas  Engineering  Section 
Mineral  Engineering  Department 
College  of  Earth  and  Mineral  Sciences 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania 

M.S.  thesis 

SUBROUTINE  ALPHAROOT 

Subroutine  uses  an  interval-halving  technique  to  find 
the  best  root  value  to  initialize  the  Newton-Raphson  (N-R) 
iterative  calculations  which  determine  the  real  root  of 
the  alpha  polynomial  on  the  interval  [0,1]. 

PARAMETERS:  delta  =  alpha  increment 

epsilon  =  alpha  convergence  criterion 
VARIABLES:  alower  =  lower  bound  of  alpha  increment 
aupper  =  upper  bound  of  alpha  increment 
f alpha  =  the  alpha  polynomial 
fprime  =  first  derivative  of  alpha  polynomial 
guess  =  iterative  variable  for  alpha 
guessO  =  initial  estimate  for  N-R 
intcount  =  #  of  intervals  until  sign  change 
iter  =  #  of  iterations  until  N-R  converged 
isign,isign2  =  flags  for  function  sign  change 
isign,isign2  =  flags  for  function  sign  change 
numroot  =  flag  for  #  of  zeros  (from  BUD  AN) 


SUBROUTINE  ALPHAROOT(jJ^comp,coefficient,xalpha,numroot, 
@  alpha) 

IMPLICIT  REAL*8(a-h,o-z) 

INTEGER  p 

DIMENSION  alpha(SOO),  coefficient(0:100),  xalpha(500) 
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PARAMETER(delta  =  0.0 ldO,  epsilon  =  l.d-06) 


* 

*  Write  table  heading 

* 


write(7,*)The  Fourier-Budan  Theorem  yields  '.numroot,'  roots  on 
@this  interval' 
write(7,3500) 

3500  format( 'Intervals ',4x,  Initial  Guess', 4x, 'Iterations', 4x, 'Calc. 

@Alpha',4x,'Exp.  Alpha 0 

Check  flag  NUMROOT  provided  by  subroutine  BUDAN  to  determine 
root-search  scheme 

if(numroot  .EQ.  0)  then 

write(6,*)  "No  root  on  the  interval  [0,1]  for  data  set  'j 
intcount  =  0 

write(7,3900)  intcount,xalpha(j) 

3900  format(i4,61x,f5.3) 

write(7,*)'No  root  on  the  interval  [0,1]' 
return 
end  if 

if(numroot  .EQ.  1)  then 
flower  =  0 
iupper  =  0 
end  if 

if(numroot  .GE.  2)  then 
ilower  =  0 
iupper  *  1 
end  if 

Use  incremental  search  to  determine  initial  guess 
Interval  Endpoint  DO-Loop 

do  0400  jroot  =  ilower,  iupper 

intcount  =  0 

Test  the  polynomial  at  endpoint  for  initial  sign  value 


if(jroot.  EQ.  ilower)  then 
guess  =  DFLOAT(ilower) 
alower  =  guess 
aupper  =  alower  +  delta 
end  if 

if(jroot  EQ.  iupper)  then 
guess  =  DFLOAT(iupper) 
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aupper  =  guess 
alower  =  aupper  -  delta 
end  if 
ichange  =  0 

0600  falpha  =  O.dO 

do  1500  p  =  1,  Ncomp 

term  =  coefficient(Ncomp-p)*guess**(Ncomp-p) 
if(  (Ncomp-p)  .EQ.  0  )  term  =  coefficient(0) 
falpha  =  falpha  +  term 

1500  continue 

Initialize  ISIGN2  on  first  pass  with  endpoint 


if(ichange  .EQ.  0)  then 
if(falpha  .GE.  0.)  then 
isign2  =  1 
else 

isign2  =  0 
end  if 
end  if 


*  Note  the  sign  of  the  function 

* 

if(falpha  .GE.  0.)  then 
isign  =  1 
else 

isign  =  0 
end  if 

*  Test  function  for  sign  change  and  increment  or  decrement  the 

*  search  variable  as  appropriate 

if(isign2  .EQ.  isign)  then 
if(jroot  .EQ.  ilower)  then 
alower  =  aupper 
aupper  =  aupper  +  delta 
guess  *  aupper 
else  if(jroot  .EQ.  iupper)  then 
aupper  =  alower 
alower  =  aupper  -  delta 
guess  =  alower 
end  if 
end  if 


♦ 

♦ 


Exit  subroutine  if  no  sign  change  is  detected  on  interval  [0,1] 
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if(  (guess  .GT.  1.)  .OR.  (guess  .LT.  0.)  )  then 
write(6,*)  'No  root  on  the  interval  [0,1]' 
write(7,3800)  intcount,xalpha(j) 
format(i4,61x,f5.3) 

write(7,*)'No  root  on  the  interval  [0,1]' 
return 
end  if 


If  NO  sign  change  but  still  within  interval,  repeat  the  sequence 

if(isign  .EQ.  isign2)  then 
isign2  =  isign 
intcount  =  intcount  +  1 
ichange  =  1 
go  to  0600 

else 

If  there  IS  a  sign  change: 

Halve  the  interval  where  the  function  crosses  the  x  axis 

guessO  =  (alower  +  aupper)  /  2.d0 
end  if 


Provide  this  guess  to  Newton-Raphson  to  begin  calculations 


guess  =  guessO 


N-R  is  limited  to  1000  iterations  for  convergence 


iter  =  0 

do  1000  iterlimit  =  1,  1000 

iter  =  iter  +  1 
f alpha  =  O.dOO 
fprime  =  O.dOO 

do  2000  p  =  1,  Ncomp 

falpha  =  falpha  +  coefficient(Ncomp-p) 
*guess**(Ncomp-p) 

fprime  =  fprime  +  (Ncomp-p)*coefficient(Ncomp~p) 
*guess**(Ncomp-p-l) 

continue 

calc  =  guess  -  falpha/fprime 
error  =  DABS((calc  -  guess)/calc) 
guess  =  calc 

if(erroi'  .LE.  epsilon)  go  to  3000 
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1000  continue 

print  *,'N-R  method  failed  to  converge  after  1000  iterations' 


Output  results  to  file  "TABLE" 


3000  write(7,3600)  intcount,guessO,iter,guess,xalpha(j) 

3600  format(j4,13x,f5.3,10x,i4,13xff9.6,7x,f5.3) 

alpha(j)  =  guess 

Begin  search  for  root  from  opposite  end  of  interval 


0400  continue 

return 

end 

************************************************************************ 

* 
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M.S.  thesis 

SUBROUTINE  ALPHAPLOT 

This  subroutine  is  used  for  several  purposes: 

1.  Plotting  F(alpha)  vs  alpha  [Rachford-Rice  obj  function] 

2.  Plotting  F(alphai  vs  alpha  [polynomial] 

3.  Plotting  Fprime  vs  alpha  [polynomial] 

*********************************************************************** 

SUBROUTINE  ALPHAPLOT(NcompJ,molefrac,alpha,coefTicient,Ki) 

IMPLICIT  REAL*8  (a-h,o-z) 

REAL* 8  Ki(500,100)jnolefrac(100) 

DIMENSION  alpha(500),coefficient(0: 100) 

INTEGER  p 

PARAMETER ( start  =  O.OdO,  end  =  2.0d0,  stepsize  =  0.0005d0) 
OPEN(unit=  1 1  ,file=  'fa.  plot  ',status=  'unknown  0 


OPEN (uni  t=  1 2,file=  'fprime. plot',status=  'unknown  0 
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* 

*  Number  of  data  points  for  plotting 


number  =  IDINT((end  —  start  +  stepsize)/stepsize) 

************************************************************************ 

*  F(alpha)  vs  alpha  [polynomial] 

*  F '(alpha)  vs  alpha  [polynomial] 

*  Adjust  NcompJVpress  in  PARAMETER 


write(ll,*)  number 
do  1000  phase  =  start, end, stepsize 
falpha  =  O.dOO 
fjprime  =  O.dOO 
do  2000  p  =  l.Ncomp 

falpha  =  falpha  +  coefficient(Ncomp-p)* 

@  phase**(Ncomp-p) 

C  fprime  -  fprime  +  ( Ncomp-p)  *  coefficient Ncomp-p )  * 

C  @  phase**  (Ncomp-p-1) 

2000  continue 

write( 11,3600)  phase,falpha 
C  write(113600)  phase  fprime 

3600  format(f7.3,2x,f25.12) 

1000  continue 

*************************************************** * * * *** *** *** * * * * ** *** 


C* 


Rachford-Rice  objective  function 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c* 

C  4000 

c 

c 

C  3500 

c* 

C  3000 


do  4500  k  =  Impress 
k  =  6 

write(lj,*)  number 

do  3000  phase  =  start, end, stepsize 
falpha  =  O.dOO 
do  4000  i  =  Iff  comp 

falpha  =  falpha  +  (molefrac(i)*(Ki(k,i)  -  l.dO)) 
@  (l.d00  +  phase*(Ki(k,i)  -  J.dO)) 

End  of  i  loop 
continue 

write(  1 1 J500)  phase  falpha 
format(f732xf25.12) 

End  af  phase  loop 

continue 


C*  End  af  k  loop 

C  4500  continue 


/ 


*********************************************************************** 


CLOSE(unit=ll) 

CLOSE(unit=12) 

return 

end 
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****************4^*******£±±***************************************** 


* 

*  6  Dec  91 

* 

*  BRETT  D.  WEIGLE 

*  Petroleum  and  Natural  Gas  Engineering  Section 

*  Mineral  Engineering  Department 

*  College  of  Earth  and  Mineral  Sciences 

*  The  Pennsylvania  State  University 

*  University  Park,  Pennsylvania 

* 

*  M.S.  thesis 

* 


*  SUBROUTINE  BUD  AN 

* 

*  Subroutine  uses  the  Fourier-Budan  Theorem  to  determine 

*  the  number  of  roots  that  the  alpha  polynomial  has  on  me 

*  interval  [u,v]. 

* 


*  PARAMETERS:  iu  =  lower  bound  of  alpha  interval 

*  iv  =  uppper  bound  of  alpha  interval 

*  VARIABLES:  coefficient  =  coefficient  of  alpha  polynomial 

*  dcoeff  =  coefficient  of  polynomial  derivatives 

*  deriv  =  derivatives  of  alpha  polynomial 

*  fvapor  =  the  alpha  polynomial 

*  ia.ib  =  #  of  sign  changes  for  derivative  series 

*  ivapor  =  alpha  =  vapor  fraction 

*  jsign.ksign  =  flags  for  derivative  sign  change 

*  numroot  =  number  of  zeros  on  the  interval 
********************************************************************* 


SUBROUTINE  BUDAN(J,Ncomp,coefficient,numroot) 

IMPLICIT  REAL*8(a-h,0-z) 

INTEGER  p 

DIMENSION  dcoeff(0:100, 0:100),  coefficient(0:100),  deriv(0:100) 
PARAMETER(iu  =  0,  iv  =  I) 


C  DATA  (cotfficientfl),  l  =  0J1comp-l)  /- l.,l.,-2.,3.,-4.J5 .  / 

OPEN(unit=2,file='test',status=  'unknown  0 
REWIND(unit=2) 


ia  =  0 
ib  =  0 

do  0500  ivapor  =  iu,  iv,  1 

Evaluate  the  polynomial  function  at  the  endpoints  iu  and  iv 


fvapor  =  O.dO 
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do  0600  p  =  1,  Ncomp 

fvapor  =  fvapor  +  coefficient(Ncomp-p)*ivapor**(Ncomp-p) 
0600  continue 

write(2,*)  'fvapor  =  '.fvapor 
write(2,*) 


Calculate  coefficients  of  first  derivative 


do  1000  n  =  Ncomp-1,  0,  -1 
dcoeff(0,n)  =  coefficient^) 
write  (2,*)  'dcoeff(0,',n,')  =  ',dcoeff(0ji) 

1000  continue 

write(2 ,*)  '  ' 

Calculate  coefficients  of  2nd-  and  higher-order  derivatives 
as  multiples  of  those  of  the  first  derivative 

do  1500  m  =  1,  Ncomp-1 

do  2000  n  =  Ncomp-m,  1,  -1 
dcoeff(m,n-l)  =  n*dcoeff(m-l,n) 
write  (2,*)  'dcoeff(',m,',',n-l,')  = 

@dcoeff(mji-l) 

2000  continue 

write( 2,*)  ' 

1500  continue 

i 

Evaluate  the  derivative  series  at  the  endpoints  iu  and  iv 


do  3000  m  =  1,  Ncomp-1 
deriv(m)  =  O.dO 

do  4000  n  =  Ncomp-m,  1,  -1 

term  =  dcoeff(m,n-l)*ivapor**(n-l) 
if(  (n-1)  .EQ.  0  )  term  =  dcoeff(m,n-l) 
deriv(m)  =  deriv(m)  +  term 
write(2,*)  'inter  deriv(',ra,')  =  ',deriv(m) 
4000  continue 

write(2,*)  'total  deriv(',ra,')  =  ',deriv(m) 
write(2,*) 

3000  continue 


Count  the  sign  changes  between  the  terms  of  the  series 


if(fvapor.  LT.  0.)  then 
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ksign  =  0 
else 

ksign  =  1 
end  if 

write(2,*)  Tcsign  =  '.ksign,'  for  fvapor' 

do  5000  i  =  1,  Ncomp-1 
if(deriv(i)  .LT.  0.)  then 
jsign  =  0 
else 

jsign  =  1 
end  if 

write(2,*)  'jsign  =  '.jsign,'  for  deriv(',i,')' 


*  Increment  A  or  B,  depending  upon  the  endpoint  under  evaluation 

* 

if(ivapor  .EQ.  iu)  then 
if(ksign  .NE.  jsign)  then 
ia  =  ia  +  1 

write(2,*)  'ia  =  '.ia,'  for  deriv(',i,y 
end  if 
end  if 

if(ivapor  .EQ.  iv)  then 
if(ksign  .NE.  jsign)  then 
ib  =  ib  +  1 

write(2,*)  'ib  =  '.ib,'  for  deriv(',i,')' 
end  if 

>  end  if 

ksign  =  jsign 

write(2,*)  Tcsign  =  '.ksign,'  after  deriv(',i,')' 
write(2,*)  '  ' 

5000  continue 
0500  continue 


* 

*  Pass  a  flag  to  calling  program  to  indicate  root  conditions 

* 


write(2,*)  'ia  =  ',ia,'  and  ib  =  '4b 

numroot  =  ia  -  ib 

write(2,*)  'numroot  =  '.numroot 

write(2,6000)  Ncomp-1,  numroot,  iu,  iv,  J 
6000  format(This  polynomial  of  order  '43,'  has  '43,'  zeros  on  the  in 
@terval  ['42, ','42,']  for  J  =  '43) 

CLOSE(unit=2) 

return 

end 


